Matlab:离散时间滤波器设计中的困惑

Matlab: Confusion in Discrete time filter design

我在matlab中有以下代码:

L = 10000;
PN_30dB = -100; % dBc per Hz at 10k
f_offset = 5e6;
f_offset2 = 10e5;
F0 = 2.5e9;
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;
sigma3 = (f_offset2/F0)*sqrt(2*10^(PN_30dB/10)/F0);
y1 = zeros(1,L);
y2 = zeros(1,L);
y3 = zeros(1,L);
y4 = zeros(1,L);
y5 = zeros(1,L);
y = zeros(1,L);
x = zeros(1,L);
for i = 2:L,
    x(i) = sigma3*randn(1);
    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
fft1 = fft(y);
fft1 = fft1(1:length(y)/2+1);
psd1 = (1/(F0*length(y)))*abs(fft1).^2;
psd1(2:end-1) = 2*psd1(2:end-1);
freq = 0:F0/length(y):F0/2;
figure(3);
semilogx(freq,10*log10(psd1))
grid on

acc = 0;
actual_timestamps_3 = zeros(1,L);
flicker = zeros(1,L);
for i = 1:L,
   acc = acc + y(i);
   actual_timestamps_3(i) = i/F0 + acc;
   flicker(i) = acc;
end

fft1 = fft(2*pi*F0*flicker);
fft1 = fft1(1:length(flicker)/2+1);
psd1 = (1/(F0*length(flicker)))*abs(fft1).^2;
psd1(2:end-1) = 2*psd1(2:end-1);
freq = 0:F0/length(flicker):F0/2;
figure(4);
semilogx(freq,10*log10(psd1))
grid on

在这段代码中,我试图创建一个名为 flicker 的输出信号,它的功率谱密度应该有 30dB/dec 的衰减。
为此,我正在累积一个名为 y(i) 的信号(在第二个 for 循环中),该信号具有 10dB/dec 的衰减,如代码中的图 3 所示。由于累积应该再增加 20db/dec,我希望闪烁信号有 30dB/dec 的滚降。
信号 y(i) 是第一个 for 循环中实现的离散时间滤波器的输出。
但是我没有看到闪烁信号的预期衰减 (30dB/dec)(图 4)。该图显示闪烁信号仅为 20dB/dec。有人可以解释我做错了什么吗? 编辑 代码中的图4如下所示:

当您使用 FFT 估计功率谱密度时,您正在查看的数据 (flicker) 有效地乘以长矩形 window。 window 引入了一些 spectral leakage,不幸的是,矩形 windows 的频谱泄漏衰减小于 20dB/decade。较强频率分量的泄漏因此掩盖了您试图观察的衰减。

为避免这种情况,您应该将信号乘以另一个 window function。有很多实验(提供不同的权衡),但出于说明目的,您可以使用 blackman window 和:

fft1 = fft(2*pi*F0*flicker .* blackman(length(flicker))');