差分方程求解、绘图及输出前四周期计算方法咨询

Hey there! Let's tackle this problem step by step— I'll cover everything from solving the difference equation, plotting the output, to calculating the first four cycles of the output directly. No overly complex jargon, just clear, actionable steps.

Your input is a 500Hz sine wave sampled at 6000Hz. Let's translate that into a discrete-time formula:

• Sampling rate ( F_s = 6000 , \text{Hz} ), input frequency ( f = 500 , \text{Hz} )
• Angular frequency per sample: ( \omega = 2\pi \frac{f}{F_s} = \frac{\pi}{6} )
• Discrete-time input: ( x(n) = \sin\left(\frac{\pi n}{6}\right) ) for ( n \geq 0 ), and ( x(n) = 0 ) for ( n < 0 )
• Each cycle of ( x(n) ) has ( \frac{F_s}{f} = 12 ) samples (so we'll need 48 total samples for 4 full cycles)

Your equation is:

( y(n) - 2.56y(n-1) + 2.22y(n-2) - 0.65y(n-3) = x(n) + x(n-3) )

Rewrite it as a recursive formula (this is the easiest way to compute outputs directly):
( y(n) = 2.56y(n-1) - 2.22y(n-2) + 0.65y(n-3) + x(n) + x(n-3) )

Key Initial Condition #

Since this is a 3rd-order equation, we need 3 initial values for ( y ). We assume the system is at rest initially, so:
( y(-1) = y(-2) = y(-3) = 0 )
And ( x(n) = 0 ) for ( n < 0 ), so ( x(n-3) = 0 ) when ( n < 3 )

Let's compute the first few values manually to show the pattern, then you can extend this to 48 samples:

• ( n=0 ): ( y(0) = 0 + 0 + 0 + x(0) + x(-3) = 0 + 0 = 0 )
• ( n=1 ): ( y(1) = 2.56y(0) - 2.22y(-1) + 0.65y(-2) + x(1) + x(-2) = 0 + 0 + 0 + \sin(\pi/6) + 0 = 0.5 )
• ( n=2 ): ( y(2) = 2.56y(1) - 2.22y(0) + 0.65y(-1) + x(2) + x(-1) = 2.56*0.5 + 0 + 0 + \sin(\pi/3) + 0 ≈ 1.28 + 0.866 = 2.146 )
• ( n=3 ): ( y(3) = 2.56y(2) - 2.22y(1) + 0.65y(0) + x(3) + x(0) ≈ 2.562.146 - 2.220.5 + 0 + 1 + 0 ≈ 5.494 - 1.11 + 1 = 5.384 )

Continue this recursive calculation up to ( n=47 ) (the end of the 4th cycle). For larger calculations, use a script (see section 4) to avoid manual errors.

Discrete signals are best visualized with stem plots. Here's a complete script to generate the input, compute the output, and plot both:

import numpy as np
import matplotlib.pyplot as plt

# Signal parameters
Fs = 6000
f = 500
omega = 2 * np.pi * f / Fs
num_cycles = 4
samples_per_cycle = int(Fs / f)  # 12 samples per cycle
total_samples = num_cycles * samples_per_cycle

# Generate input signal x(n)
n = np.arange(total_samples)
x = np.sin(omega * n)
# Create delayed input x(n-3) (pad first 3 samples with 0)
x_delay3 = np.concatenate([np.zeros(3), x[:-3]]) if total_samples >3 else np.zeros(total_samples)

# Initialize output array
y = np.zeros(total_samples)

# Recursively compute y(n)
for i in range(total_samples):
    # Handle initial conditions (use 0 for out-of-bounds indices)
    y_prev1 = y[i-1] if i >=1 else 0
    y_prev2 = y[i-2] if i >=2 else 0
    y_prev3 = y[i-3] if i >=3 else 0
    x_delay = x_delay3[i] if i >=3 else 0
    
    y[i] = 2.56*y_prev1 - 2.22*y_prev2 + 0.65*y_prev3 + x[i] + x_delay

# Plot results
plt.figure(figsize=(12, 7))

# Input signal plot
plt.subplot(2, 1, 1)
plt.stem(n, x, basefmt="b-", linefmt="b-", markerfmt="bo")
plt.title("Input Signal x(n) (500Hz, 6000Hz Sampling)")
plt.xlabel("Sample Index n")
plt.ylabel("Amplitude")
plt.xlim(-0.5, total_samples - 0.5)

# Output signal plot
plt.subplot(2, 1, 2)
plt.stem(n, y, basefmt="r-", linefmt="r-", markerfmt="ro")
plt.title("Output Signal y(n) from Difference Equation")
plt.xlabel("Sample Index n")
plt.ylabel("Amplitude")
plt.xlim(-0.5, total_samples - 0.5)

plt.tight_layout()
plt.show()

What This Script Does: #

1. Generates the input sine wave and its 3-sample delayed version
2. Uses a loop to compute each ( y(n) ) with proper handling of initial conditions
3. Plots both signals with stem plots (ideal for discrete-time data)

If you want a mathematical formula for ( y(n) ), use the Z-transform:

1. Take the Z-transform of both sides of the difference equation
2. Derive the system function ( H(z) = \frac{Y(z)}{X(z)} = \frac{z^3 + 1}{z^3 - 2.56z^2 + 2.22z - 0.65} )
3. Multiply by the Z-transform of ( x(n) ), then perform partial fraction expansion to find the inverse Z-transform. This is more complex, but useful for understanding system behavior.

内容的提问来源于stack exchange，提问作者Zuzu