2017-07-30 101 views
2

我在MATLAB以下代码:混乱在MATLAB滤波器滚降

L = 10000; 
Power = -100; 
A = 10^0.5; 
Fs = 25e6; 
fc1 = 100; 
fc2 = 1e3; 
fc3 = 10e3; 
fc4 = 100e3; 
fc5 = 1e6; 
a1 = 2*pi*fc1/Fs; 
a2 = 2*pi*fc2/Fs; 
a3 = 2*pi*fc3/Fs; 
a4 = 2*pi*fc4/Fs; 
a5 = 2*pi*fc5/Fs; 

x = wgn(1,L,Power); 
y1 = zeros(1,L); 
y2 = zeros(1,L); 
y3 = zeros(1,L); 
y4 = zeros(1,L); 
y5 = zeros(1,L); 
y = zeros(1,L); 
for i = 2:L, 
    y1(i) = (1-a1)*y1(i-1) + a1*x(i); 
    y2(i) = (1-a2)*y2(i-1) + a2*x(i)/A; 
    y3(i) = (1-a2)*y3(i-1) + a3*x(i)/A^2; 
    y4(i) = (1-a2)*y4(i-1) + a4*x(i)/A^3; 
    y5(i) = (1-a2)*y5(i-1) + a5*x(i)/A^4; 
    y(i) = y1(i) + y2(i) + y3(i) + y4(i) + y5(i); 
end 
fft1 = fft(y); 
fft1 = fft1(1:length(y)/2+1); 
psd1 = (1/(Fs*length(y)))*abs(fft1).^2; 
psd1(2:end-1) = 2*psd1(2:end-1); 
freq = 0:Fs/length(y):Fs/2; 
figure(3); 
semilogx(freq,10*log10(psd1)) 
grid on 

Ts = 40e-9; 
z = tf('z',Ts); 
H1 = a1/(1-(1-a1)*z^-1); 
H2 = (a2/A)/(1-(1-a2)*z^-1); 
H3 = (a3/A^2)/(1-(1-a3)*z^-1); 
H4 = (a4/A^3)/(1-(1-a4)*z^-1); 
H5 = (a5/A^4)/(1-(1-a5)*z^-1); 
H = (H1 + H2 + H3 + H4 + H5); 
figure(5); 
bode(H),grid 

该代码的意图是闪烁噪声,它有在它的功率谱密度斜率10dB的/癸建模。
若要建模,有一个过滤器的输出是y(i)和输入x(i)这是白色高斯噪声在这种情况下。在滤波器的波特图中(标记为图5),我可以看到它具有10dB/dec的滚降速度。
但是,当我检查输出噪声(在这种情况下y(i))功率谱密度(标记为图3)时,我看到20dB/dec的滚降。
有人能解释我做错了什么,为什么我无法获得10dB的功率谱密度滚降?

+0

它看起来像这将有助于让你通过调试运行此。 – beaker

回答

1

您的波特图使用的传递函数是5个相似的子滤波器的总和,每个滤波器都使用系数a i。在你的实现中,你应该同样有5个具有规则结构的子过滤器。 Howerver在复制每一行时,已将y3,y4y5的一些系数保留为使用a2。你应该简单地替代了相应的a3a4a5像这样让您预期的结果:

for i = 2:L, 
    y1(i) = (1-a1)*y1(i-1) + a1*x(i); 
    y2(i) = (1-a2)*y2(i-1) + a2*x(i)/A; 
    y3(i) = (1-a3)*y3(i-1) + a3*x(i)/A^2; 
    y4(i) = (1-a4)*y4(i-1) + a4*x(i)/A^3; 
    y5(i) = (1-a5)*y5(i-1) + a5*x(i)/A^4; 
    y(i) = y1(i) + y2(i) + y3(i) + y4(i) + y5(i); 
end