Zero Padding Tutorial M-file

This is the code I used in the youtube tutorial on zero padding (see http://youtu.be/7yYJZfc5q-I).

% This tutorial attempts to explain how zero-padding works.
% You should already have some exposure to how the DFT works and
% be have seen the effects of zero paddding and windowing.
%

fs = 1000;
T = 1/fs;
N = 1000; % Create a signal with N samples = 1 second
n = 0:N-1; % sample numbers
t = n*T; % sampling times over a one second duration.
x = cos(2*pi*3*t);

%plot the magnitude spectrum 
X_mags = abs(fft(x));

subplot(2,1,1)
plot(x)
xlabel('Samples');
ylabel('Amplitude')
title('Time-Domain Signal');

num_disp_bins = 15;
subplot(2,1,2)
plot([0:num_disp_bins-1], X_mags(1:num_disp_bins));
hold on
plot([0:num_disp_bins-1], X_mags(1:num_disp_bins),'k.');
hold off
xlabel('Frequency Bins');
ylabel('Magnitude');
title('Frequency Content Magnitudes');

%Illustrate Basis functions
plot(n,x,'k')
xlabel('Samples');
ylabel('Amplitude')
for k = 0 : 20 
    plot(x,'kx')
    hold on
    cos_basis = cos(2*pi*k*n/N);
    sin_basis = sin(2*pi*k*n/N);
    
    plot(n,cos_basis,'r')
    plot(n,sin_basis,'g')
    title({['Signal being anlaysed (black) with basis functions'...
        ' that have ' num2str(k) ' cyles over '  ...
        num2str(N) ' samples']...
        ['Cosine basis function shown in red, Sine in green']})
    hold off
    pause
end

Advertisement

PDF Version of “Frequency Analysis using the DFT – A practical example”

PDF Version of “Frequency Analysis using the DFT – A practical example”

This PDF was generated using matlab’s publishing tool on the code available at https://dadorran.wordpress.com/2013/12/17/frequency-analysis-using-the-dft-a-practical-example/

Frequency Analysis using the DFT – A practical example

%% Frequency Analysis using the DFT - A practical example
% This tutorial shows how to analyse the frequency content of a
% bass guitar to determine the fundamental frequencies of the 
% notes being played.
% Youtube Video available at http://youtu.be/eg8eebQPfAo95
%
% Matlab code will be available at https://dadorran.wordpress.com

% First download http://eleceng.dit.ie/dorran/matlab/bass.wav
[ip fs] = wavread('bass.wav'); 
%%
% Now extract a short segment from the guitar signal 
% which contains two notes being played over 0.45 seconds
seg = ip(670000:689000);
sound(seg, fs); %use this to hear the two notes being played

plot(seg);
xlabel('samples');
ylabel('Amplitude');
%%
% Annotate the figure to show where the two notes occur and 
% indicate fundamental frequencies 
annotation('doublearrow',[0.14 .45],[0.3 0.3]);
annotation('textbox',[0.14 .25 0.4 0.05],'string',...
    {'Note 1' 'Fundamental approx 72 Hz'},...
    'LineStyle', 'none', 'HorizontalAlignment','center');
annotation('doublearrow',[0.55 .85],[0.3 0.3]);
annotation('textbox',[0.55 .25 0.4 0.05],'string',...
    {'Note 2' 'Fundamental approx 86 Hz'},...
    'LineStyle', 'none', 'HorizontalAlignment','center');

%% Extract a frame/window of length 2000 samples from Note 1   
% Analyse the frequency content of the frame.
N = 2000;
frame = seg(1:N-1); % this is a frame associated with Note 1
ft_mags = abs(fft(frame));
num_bins_to_display = 250;
figure
plot([0:num_bins_to_display-1], ft_mags(1:num_bins_to_display))
xlabel('Frequency Bins')
ylabel('Magnitude')
%%
% This plot shows that the signal has three 'strong' frequency
% componentsi.e. the fundamental and first two harmonics
%%
% Let's try to determine the fundamental frequency by analysing 
% the DFT.
% The fundamental is the strongest frequency component present 
% i.e. the maximum. By zooming in on the plot it can be seen that 
% it is located at bin number 3. Alternatively we can use matlabs 
% built-in max function.
[max_val fundamental_location] = max(ft_mags); 
% NOTE : fundamental_location uses matlabs indexing which is 
% different to bin number by 1
%%
% Each bin is separated by fs/N Hz so the fundamental frequency 
% determined from analysis of this magnitude spectrum is 3(fs/N)
% NOTE : the fundamental_location variable obtained from the max 
% function uses matlabs indexing which is different to the bin 
% number by 1.
fundamental_frequency = (fundamental_location-1)*fs/N 
%%
% This result of 66.15 Hz is inaccurate (the fundamental 
% frequency is actually about 72 Hz). 
% The inaccuracy arises due to the limitation of the DFT frequency
% resolution. 
% We can only be sure that the result obtained using this type of 
% analysis is accurate to within fs/N = 22.05 Hz.

%% Extract a frame/window of length 5000 samples from Note 1
% This increases frequency resolution and improves the accuracy of
% the analysis.
N = 5000;
frame = seg(1:N-1);
ft_mags = abs(fft(frame));
num_bins_to_display = 250;
plot([0:num_bins_to_display-1], ft_mags(1:num_bins_to_display))
xlabel('Frequency Bins')
ylabel('Magnitude')
title('Magnitude of First 250 DFT bins of Note 1')
%%
% Notice that the fundamental and harmonics are separated by a 
% greater number of DFT bins than for N = 2000.
%%
% Now, determine fundamental frequency using the same technique 
% as before.
[max_val fundamental_location] = max(ft_mags);
fundamental_frequency = (fundamental_location-1)*fs/N
%%
% The result of 70.56 Hz is more accurate because the DFT 
% frequency resolution is higher. For this case we can be sure 
% that the result will be within fs/N = 8.82 Hz of the actual 
% fundamental frequency.

%% Extract a frame/window of length 5000 samples from Note 2
N = 5000;
%extract  samples from the end of the segment
frame = seg(end-N:end); 
ft_mags = abs(fft(frame));
num_bins_to_display = 250;
figure;
plot([0:num_bins_to_display-1], ft_mags(1:num_bins_to_display))
xlabel('Frequency Bins')
ylabel('Magnitude')
title('Magnitude of First 250  DFT bins of Note 2')
%%
% Determine fundamental frequency using the same technique as above
[max_val fundamental_location] = max(ft_mags);
fundamental_frequency_note2 = (fundamental_location-1)*fs/N
%%
% This result is close to the actual fundamental frequency of 86Hz. 
% Further improvement will be obtained if N is increased further 
% (as shown in the next section)

%% Analyse the DFT of the entire segment
% The resolution is higher than previous analysis
N = length(seg);
ft_mags = abs(fft(seg));
num_bins_to_display = 250;
figure;
plot([0:num_bins_to_display-1], ft_mags(1:num_bins_to_display))
xlabel('Frequency Bins')
ylabel('Magnitude')
title('Magnitude of First 250 DFT bins of both notes')
%%
% This code might be tricky to follow - I'm just highlighting
% the fundamental and first two harmonics of each note in the plot.
% don't worry if this part of the code is confusing.
note1_fund_freq = 72;
note2_fund_freq = 86;
ft1 = ones(1, length(ft_mags))*NaN;
ft2 = ones(1, length(ft_mags))*NaN;
for k = 1: 3 
   ft1(round(note1_fund_freq*k/fs*N)...
       + 1 - 3:round(note1_fund_freq*k/fs*N)+ 1+ 3)...
       = ft_mags(round(note1_fund_freq*k/fs*N)+1 - 3:...
       round(note1_fund_freq*k/fs*N)+ 1+ 3);
   ft2(round(note2_fund_freq*k/fs*N)...
       + 1 - 3:round(note2_fund_freq*k/fs*N)+ 1+ 3)...
       = ft_mags(round(note2_fund_freq*k/fs*N)+1 - 3:...
       round(note2_fund_freq*k/fs*N)+ 1+ 3);
end
hold on
plot([0:num_bins_to_display-1], ft1(1:num_bins_to_display),'r')
plot([0:num_bins_to_display-1], ft2(1:num_bins_to_display),'g')
legend('','Note 1 fundamental and harmonics',...
    'Note 2 fundamental and harmonics')
%%
%Determine fundamental frequency using the same technique
% as above.
[max_val fundamental_location] = max(ft_mags);
fundamental_frequency_note1 = (fundamental_location-1)*fs/N
%%
% From the plot it can be seen that the frequency of the 
% fundamental of the second note is associated with bin 37 
% (no need to subtract 1 in equation below because the max 
% function isn't being used).
fundamental_frequency_note2 = 37*fs/N

%% A 'Better' Approach to higher frequency resolution 
% 
% While improvements in results have been obtained using a longer 
% window the plot of the magnitude spectrum is becoming more 
% difficult to interpret as there is a lot of spectral energy in 
% the signal (energy from two notes rather than one). It would be 
% nicer to be able to get higher frequency resolution without the 
% problem of the magnitude spectrum becoming more difficult to 
% interpret.
%
% The way to get improved resolution without having a more 
% complicated spectrum to deal with is to 'artificially' make 
% the signal longer through a process known as zero padding 
% (see https://www.youtube.com/watch?v=V_dxWuWw8yM for a tutorial
% on zero padding and windowing).
%
% We'll go back to the case where we extracted the first 5000 
% samples of note 1 and then append 40,000 samples of zero 
% amplitude to this signal.
N = 5000;
frame = seg(1:N-1);
% Make the frame 'longer' by zero-padding 40000 samples of 
% amplitude zero
zpad_frame = [frame ;zeros(40000,1)]; 
new_frame_len = length(zpad_frame);
%%
figure
plot(zpad_frame)
xlabel('samples')
ylabel('Amplitude')
title(...
 'First 5000 samples of Note 1 zero-padded by 40000 zero samples')
%%
% Now plot the magnitude spectrum as before
ft_mags = abs(fft(zpad_frame));
num_bins_to_display = 250;
figure
plot([0:num_bins_to_display-1], ft_mags(1:num_bins_to_display))
xlabel('Frequency Bins')
ylabel('Magnitude')
title('Magnitude of First 250 DFT bins of Note 1')
%%
% Now, determine fundamental frequency using the same technique 
% as before.
% NOTICE that the frequency resolution is now fs/new_frame_len
[max_val fundamental_location] = max(ft_mags);
fundamental_frequency = (fundamental_location-1)*fs/new_frame_len
%%
% This result is accurate to fs/new_frame_len Hz
fs/new_frame_len
%%
% If the frame was zero padded by a greater amount then the 
% frequency resolution would improve and therefore the accuracy of 
% the frequency estimate would improve.

%% What happens if the frame length is too short??
% We've seen in the examples so far that the separation between 
% fundamental and harmonics increases as the frequency resolution 
% increases. It follows that the separation between harmonics will
% decrease as the the frequency resolution decreases. If the 
% frequency resolution can be so low that the harmonics will not 
% be separated.
%
% Let's see what happens if we take frame length of 1500 samples
N = 1500;
frame = seg(1:N-1); % this is a frame associated with Note 1
ft_mags = abs(fft(frame));
num_bins_to_display = 100;
figure
plot([0:num_bins_to_display-1], ft_mags(1:num_bins_to_display))

xlabel('Frequency Bins')
ylabel('Magnitude')
title('Magnitude of First 100 DFT bins of Note 1')
%%
% Notice that only the first 50 bins are shown in this case 
% - this is to make it easier to see the harmonics.
%
% It can be seen that the harmonics are very close together 
% - in fact they are only separated by 1 bin!
%
% If a lower frequency resolution analysis is performed the 
% harmonics will not be separated at all.
%%
%
% Let's see what happens if we take frame length of 1200 samples
N = 1200;
frame = seg(1:N-1); % this is a frame associated with Note 1
ft_mags = abs(fft(frame));
num_bins_to_display = 100;
figure
plot([0:num_bins_to_display-1], ft_mags(1:num_bins_to_display))

xlabel('Frequency Bins')
ylabel('Magnitude')
title('Magnitude of First 100 DFT bins of Note 1')
%%
% At this stage there is now separation between the higher 
% harmonics
%
% Let's see what happens if we take frame length of 1000 samples
N = 1000;
frame = seg(1:N-1); % this is a frame associated with Note 1
ft_mags = abs(fft(frame));
num_bins_to_display = 100;
figure
plot([0:num_bins_to_display-1], ft_mags(1:num_bins_to_display))

xlabel('Frequency Bins')
ylabel('Magnitude')
title('Magnitude of First 100 DFT bins of Note 1')
%%
% There is no separation between fundamental or harmonics 
% - the frequency resolution is too low. Looking at this magnitude
% spectrum you would probably interpret it to mean that there was 
% one sinusoidal component present in the signal (or perhaps two 
% given the small peak at bin 10)
%%
% To the uninitiated the solution to this problem might be to zero
% pad the frame in order to improve the frequency resolution.
% However, this doesn't work. An indepth discussion as to reason 
% why it doesn't work is beyond the scope of this tutorial (as it
% requires a very good understanding of how the DFT works) but a 
% basic conceptual  understanding can be obtained by considering 
% an extreme case whereby you extracted a frame of which is just 
% of length 1 sample. If you only have only sample it would be 
% impossible to determine any frequency domain information since 
% frequency relates to the rate of change of samples. Even if you 
% had two samples you still have a very limited view of the 
% frequency of the signal and no amount of zero padding will help. 
% Hopefully you can appreciate that you need some number of 
% time-domain samples in order to extract out useful frequency 
% content. As a general rule you need at least one cycle of a 
% periodic signal to in order to be able to extract 'useful' 
% frequency content associated with the signal. It is only when 
% you have a sufficient amount of time-domain data that 
% zero-padding can be used to obtain more accurate frequency 
% estimates.