How to plot and normalize a spectrum in matlab

Question

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:

Answer 1

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.

Answer 2

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.