Matlab: Confusion in Discrete time filter design

Question

I have the following code in 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

In this code, I am trying to create an output signal called flicker which should have a 30dB/dec roll-off of its power spectral density.
For that I am accumulating a signal called y(i) (in the second for loop) which has a 10dB/dec roll-off as can be seen in figure 3 inside the code. Since accumulation should add another 20db/dec, I am expecting the flicker signal to have a 30dB/dec roll off.
Signal y(i) is the output of a discrete time filter implemented in first for loop.
But I am not seeing the expected roll-off (30dB/dec) for the flicker signal (in figure 4). The plot shows only 20dB/dec for the flicker signal. Could someone please explain what I am doing wrong? EDIT Figure 4 in the code is shown below:

SleuthEye · Accepted Answer · 2017-08-01T12:11:31.080

2

As you estimate the power spectral density with the FFT, you are looking at the data (flicker) which is effectively multiplied by a long rectangular window. This window introduces some spectral leakage, and unfortunately the decay of this spectral leakage for rectangular windows is less than 20dB/decade. The leakage from the stronger frequency components is thus masking the decay you are trying to observe.

To avoid this, you should multiply your signal by a different window function. There are plenty to experiment with (which offer different tradeoffs), but for illustration purposes you could use a blackman window with:

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

edited Aug 01 '17 at 12:11

answered Aug 01 '17 at 07:17

SleuthEye

14,379
2
32
61

Thank you @SleuthEye for your answer...But I was asking about the signal `flicker` which is showing 20dB/dec instead of 30dB/dec roll off. Could you please explain that. – sarthak Aug 01 '17 at 08:02
Same applies to `flicker`. When you try to get its FFT with an effective rectangular window, the leakage with less than 20dB/decade masks the 30dB/decade decay of `flicker`. – SleuthEye Aug 01 '17 at 10:30
could you please help with this question: https://stackoverflow.com/questions/46501409/power-spectral-density-calculation-in-matlab. Thanks! – sarthak Sep 30 '17 at 15:16

Matlab: Confusion in Discrete time filter design

1 Answers1