Sunday, May 9, 2021

Why is it a bad idea to filter by zeroing out FFT bins? (2020)

This question has also confused me for a long time. @hotpaw2's explanation is good. You may be interested in the simple experiment using matlab.

https://poweidsplearningpath.blogspot.com/2019/04/dftidft.html


updated information.

To verify this fact is simple, we just need to cautiously observe the spectrum of impulse response of an ideal(?) band pass filter which just zeros out FFT bins. Why do I need to add the adverb "cautiously"? If we just use the same size of the FFT to observe the response of the impulse, we will be deceived as shown in Fig 1. Nonetheless, if we add the order of DFT when observing the output of the filter, that is, zero padding the impulse response, we can find the so called Gibbs phenomenon, ripples in frequency domain, as depicted in Fig 2.

The results in fact comes from the windowing effect. If you want to entirely understand the problem, please refer to chapter 7.6 and chapter 10.1-10.2 of the bible of DSP (1). To sum up, three key points are noted here.

  1. Size of window and order of DFT(FFT) are totally independent. Don't mix them together.
  2. Properties of window (type/size) dominate the shape of DTFT. (ex. wider main lobe lead to wider transient band in frequency response.)
  3. DFT is just the sampling of DTFT in frequency domain. Moreover, the higher the order of DFT, the denser the spectrum of DFT is.

So, with the help of denser spectrum in Fig.2, we can see through the mask of ideal(fake) Band pass filter.

enter image description here Deceitfully Freq. Response.

enter image description here Gibbs phenomenon in Freq. Response.

(1) Alan V. Oppenheim and Ronald W. Schafer. 2009. Discrete-Time Signal Processing (3rd ed.). Prentice Hall Press, Upper Saddle River, NJ, USA.

fps = 15;

LPF = 1;
HPF = 2;

n = -511:512;
n0 = 0;
imp = (n==n0);

NyquistF = 1/2*fps;

%% Ideal BPF
tmp_N = 512;
tmp_n = 0:1:tmp_N-1;
freq = ( n .* fps) ./ tmp_N;
F = fft(imp, tmp_N);  
F_bpf = IdealBandpassFilter(F, fps, LPF, HPF);
imp_rep =[real(ifft(F_bpf))'];

% Zero padding.
imp_rep2 =[zeros(1,2048) real(ifft(F_bpf))' zeros(1,2048)];

N = 2^nextpow2(length(imp_rep));
F = fft(imp_rep,N);
freq_step = fps/N;
freq = -fps/2:freq_step:fps/2-freq_step;
freq = freq(N/2+1:end)';

figure;
plot(freq,abs(F(1:N/2)));
xlabel('freq(Hz)');
ylabel('mag');
title('Mis leading Freq Response');


N = 2^nextpow2(length(imp_rep2));
F = fft(imp_rep2,N);
freq_step = fps/N;
freq = -fps/2:freq_step:fps/2-freq_step;
freq = freq(N/2+1:end)';

figure;
plot(freq,abs(F(1:N/2)));
xlabel('freq(Hz)');
ylabel('mag');
title('Zero Padding (DFT) with more points');

%% Function
function filered_signal = IdealBandpassFilter(input_signal, fs, w1, w2)

    N = length(input_signal);
    n = 0:1:N-1;
    freq = ( n .* fs) ./ N;

    filered_signal = zeros(N, 1);

    for i = 1:N
        if freq(i) > w1 & freq(i) < w2
            filered_signal(i) = input_signal(i);
        end

    end
end


from Hacker News https://ift.tt/3ey6P9G

No comments:

Post a Comment

Note: Only a member of this blog may post a comment.