1
$\begingroup$

I am given the coefficients of an FIR filter. I have created its impulse response through a convolution with a unit pulse. Now I want to plot the frequency spectrum. I perform zero-padding and windowing with a blackman-window. This does not seem to have a big effect on my result though. Still, I wonder: For some reason my amplitude of the resulting spectrum is not normalized correctly. Can somebody see what I am doing wrong? 

fs = 44100;%sampling frequency
Ts = 1/fs; %time step
t = 0:1/fs:2-Ts; %signal duration

unit_pulse = zeros(length(t),1);
unit_pulse(1) = 1;

b = [1/8 1/8 1/8 1/8 1/8 1/8 1/8 1/8];
ir_filter = filter(b,1,unit_pulse);

nfft = length(ir_filter); %length of the time domain signal
%in order to obtain a good frequency resolution, pad the signal with zeros
%such that the length of the signal is a power of 2
no_zeros = (2^nextpow2(nfft)-nfft)/2;
padding = zeros(no_zeros,1);
signal = vertcat(padding, ir_filter, padding);
windowed_signal = signal .* blackman(length(signal));

spectrum_2side = real(fftshift(fft(windowed_signal))) ./ nfft;
f_2side = linspace(-fs/2, fs/2, length(windowed_signal));
spectrum_1side = 2*spectrum_2side(end/2+1:end);
f_1side = f_2side(end/2+1:end);

fig2 = figure();
hold on
grid on
plot(f_1side, spectrum_1side, 'b');
plot(f_2side, spectrum_2side, 'r');

This is a figure of what the impulse response looks like: stemplot of ir_filter

Improve this question
$\endgroup$
0
0
$\begingroup$

Luk the following modified line corrects the DFT magnitude scale problem:

 signal = vertcat(padding, fftshift(ir_filter), padding);

While concatenating the filter impulse response into the zero paded signal, you loose the symmetry and the consequent multiplication with the symmetrical window (Blackman forexample) therefore reduces your filter much more than it should. By using the fftshift there, I'm guiding the zero padded signal into a symmetric final form.

Note also that, when you generate the impulse response via the line 

ir_filter = filter(b,1,unit_pulse);

it won't give you a symmetrical result.

Finally, your spectrum computation is also wrong and the following is the corrected line: 

 spectrum_2side = abs(fftshift(fft(windowed_signal)));

If you want to see the conventional DTFT / DFT reslt, then do not divide by the DFT length. Also use the magnitude for display, unless you want to inspect th real part exclusively.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Hey, thx very much for your help! I ran the code by changing signal into vsignal. But my amplitude of the frequency response is still not = 1, which I expect it to be. You are saying that my impulse response is not symmetric. And you're right about that, it is just 8 pulses with all of them on the positive x-axis. Why is it important to have a symmetric signal, and how could I create one? $\endgroup$ – user38624 Nov 5 '18 at 12:43
  • $\begingroup$ Hi! vsignal was a typo :-) you don't have to change the name of the signal. Center of your filter and center of the window should match for proper windowing. Finally, your spectrum computation is also wrong. And I'll show you the correction in the answer. $\endgroup$ – Fat32 Nov 5 '18 at 13:00
  • $\begingroup$ cool, now it's working! Thx so much! ... still: Can you explain to me why your modification of my signal has made it symmetric? ir_filter has 88200 entries, only the first 8 are non-zero and have a value of 0.125. To my knowledge, fftshift() stacks entries (end/2:end) onto the beginning of the vector. So fftshift(ir_filter) has 44100 zeros in the beginning of the vector, then 8 entries of 0.125 and 44092 zeros at the end of the vector. Ok, you're right, it is almost symmetric. $\endgroup$ – user38624 Nov 5 '18 at 20:18
  • $\begingroup$ And you are saying that this is important, because the blackman window is symmetric? Why do you nor normalize the spectrum? Shouldn't it be normalized? Otherwise the amplitude of the spectrum will become greater the more non-zero entries my signal has. $\endgroup$ – user38624 Nov 5 '18 at 20:19
  • $\begingroup$ You are right @Luk ! you know this thing better than me ;-) It's not symmetric it's just almost symmetric. I have lazily abused fftshift() to make your impulse response almost symetric. Since the length was soo long, the inaccuracy is small. To learn more about the topic please read FIR filter design window method ... $\endgroup$ – Fat32 Nov 5 '18 at 20:30
0
$\begingroup$

Short answer::

Use 

freqz(coeff)

Long answer:

To plot the frequency response for a filter (which is a continuous function) you would use the Discrete Time Fourier Transform (DTFT), not the Discrete Fourier Transform (DFT) which is a discrete function. The FFT is an algorithm that implements the DFT. The difference between the DTFT and the DFT is that the DTFT is a summation from $-\infty% to $+\infty% while the DFT is a summation only over N samples. This is given in the formulas below.

DTFT: $$X(e^{j\omega}) = \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n}$$

DFT: $$X(k) = \sum_{n=0}^{N-1}x[n]e^{-j\omega_o n}$$

with k = 0,1,2 ... N-1 and $\omega_o = \frac{2\pi}{N}$

The DFT is a sampled version of the DTFT (the samples of the DFT are on the DTFT). You can approximate the DTFT using the DFT with zero padding (as that emulates extending the time domain to $\pm \infty$); the more zeros you add, the more samples will be interpolated in the frequency domain, and all samples will be samples on the DTFT. (So not just the next power of 2 as the OP has done, which is done to use the fft algorithm more efficiently, but enough to visually interpolate the entire spectrum. Note this does NOT increase frequency resolution but plots more samples on the underlying DTFT).

The freqz in Matlab/Octave and Python command also given above returns the frequency response and is the simplest approach.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy