[Ch7. Filter Design] Why is it a bad idea to filter by zeroing out FFT bins? (為何不能DFT轉過去,直接設成想要的形狀,再IDFT轉回來?)

The complexity of DFT/IDFT using FFT/IFFT is only n*log2n whereas the complexity of convolution is n square. Then, why do we need to learn so many methods to design a filter? Why is it a bad idea to filter by zeroing out FFT bins? This question might be the most confusing one when learning DSP.

To verify this fact is simple, we just need to cautiously observe the spectrum of the impulse response of an ideal(?) bandpass 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 the frequency domain, as depicted in Fig.2.

The results in fact come 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 leads to the wider transient band in frequency response.)
3. DFT is just the sampling of DTFT in the frequency domain. Moreover, the higher the order of DFT, the denser the spectrum of DFT is.

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


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



這大概是DSP最容易令初學者困惑的問題之一。
乍看之下有道理,但因為window effect 導致事與願違。
其實只要使用利用DFT是DTFT於頻域的sample此特性,多去sample一下(zero padding),即可看出端倪。
(可參考聖經第七章FIR with Window & 10.1-10.2)

Fig.1. Deceitfully Freq. Response.

Fig.2. Gibbs phenomenon in Freq. Response.
.
-
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


Comments

Popular posts from this blog

[Ch4 ADC/DAC] How to simulate ADC/DAC process in Matlab? How to actually reconstruct a signal nearly sampled in Nyquist Rate?

[Chapter 5.1] Meaning of General Linear Phase or Group Delay of a filter?