Plotting Data

2D Plots

The following table summarizes the list of functions provided by MATLAB for the creation and customization of two-dimensional plots.

MATLAB Plot Functions

Function	Description
plot	creates a plot
loglog	creates a plot with log scales for both axes
semilogx	uses a log scale for the x-axis
semilogy	uses a log scale for the y-axis
title	adds a title to the plot
xlabel	adda a label to the x-axis
ylabel	adds a label to the y-axis
text	displays a text string at a specified location
gtext	same as text but allows use of the mouse to position text
grid	turns on grid lines

The plot function can accept one, two, or more arguments and produces a plot of the data contained in the arguments. When a single vector argument is passed to plot, the elements of the vector form the dependent data and the index of the elements form the dependent data. For example the following sequence of commands:

>> x=[1 4 9 16 25 36 49 64 81 100]
>> plot(x)

results in the plot displayed in Figure 1. The output is a very simple graph that lets you painlessly view your data. A specific set of values for the independent axis can be passed to the plot function via additional arguments. It is also possible to make multiple plots on a single graph. For example, if a plot of the functions x=sin(t) and y=cos(t) for values of t ranging from -p to +p, the following sequence of commands could be used:

>> t=-pi:pi/100:pi;
>> x=sin(t);
>> y=cos(t);
>> plot(t,x,'r-',t,y,'g--');
>> title('t vs. sin(t) and cos(t)')
>> xlabel('t')
>> ylabel('sin(t) and cos(t)')

[Missing image] Figure 1. Simple Plot

The resulting plot is shown in Figure 2.
[Missing image]
Figure 2. Multiple Plots

There are a some points of interest in the sequence of commands displayed above:

• The semi-colons at the end of some of the commands above indicate to MATLAB that the result of the command should not be echoed to the screen.
• The first command generates a vector of values between the range [-p,p] increments of p/100 and assigns the name t to the vector
• The plot command makes plots of t versus sin(t) using a solid red line (r-) and t versus cos(t) using a dashed green line (g--).
• The title, xlabel, and ylabel commands simply specify the title, x-axis label and y-axis labels of the graph.

The arguments to plot enclosed in quotes specify the type and color of the line used to draw the line connecting the datapoints. Table 3 below lists the colors and symbols that are provided by MATLAB:

MATLAB Plot Symbols

Symbol	Color	Symbol	Linestyle
y	yellow	.	point
m	magenta	o	circle
c	cyan	x	x-mark
r	red	+	plus
g	green	*	star
b	blue	-	solid
w	white	:	dotted
k	black	-.	dash-dot
		--	dashed

The plot function can also accept matrices as arguments. The are four cases that the user should be aware of when passing matrices as arguments to the plot function. The four cases and their ramifications are:

2. A single matrix argument is passed to plot. In this situation, each column of the matrix is treated as a dependent dataset and the row number is treated as the independent dataset.

3. Two arguments with the first argument being a vector, x, and the second argument being a matrix ,Y , are passed to plot. The rows OR columns of Y are plotted versus the vector x. The selection of wether the rows or columns should be used depends on the size of the vector x. If Y is a square matrix, then its columns are used.

4. Same as in case 2 above except that the order of the arguments are reversed. In this case the vector x is plotted versus each row or column of Y.

5. If a matrix, X, and a matrix, Y, are passed as arguments to the plot function, then the columns of Y are plotted versus the columns of X.

The following examples illustrate each of the cases listed above.
Example 1: X is a matrix argument to plot:
[Missing image]

>>X=[1 2 3; 4 5 6; 7 8 9];
>>plot(X);

Figure 3. Example 1







Example 2: x is a vector and is the first argument to[Missing image] plot, and Y is a matrix and is the second argument to plot:

>>x=[2 4 6 8 10];
>>Y=[3 5 7 9 11; 4 8 12 16 20];
>>plot(x,Y,'*');

Figure 4. Example 2
[Missing image]Example 3: x is a vector and is the second argument to plot, and Y is a matrix and is the first argument to plot - the data from Example 2 is used in this example.:

>>plot(Y,x,'o');

Figure 5. Example 3




Example 4: both X and Y are matrices and are arguments to plot.[Missing image]

>>X=[1 2 3; 4 5 6; 7 8 9];
>>Y=[2 4 6; 16 25 36; 10 11 12];
>>plot(X,Y);

Figure 6. Example 4

Plotting Data Contained in a Data File

Frequently, you will need to plot data that has been saved in a file. This data may have been generated by a program or entered by a user by hand. In such an instance, you can use the load command to read in the data from the file.As an example, consider a data file called rungKut.dat that at contains three columns of data with the first column representing the independent dataset and the remaining two columns representing 2 dependent datasets. The first step is to read the data into a matrix. Once the data has been loaded into a matrix, plotting of the data can be performed in the same manner as was discussed in the previous section. The load command automatically assigns the data contained in file to a matrix with the same name as the file. The following example illustrates how the data from rungKat.dat is read into MATLAB and the graphed:

>> load rungKut.dat
>> plot(rungKut(:,1), rungKut(:,2), 'g-', rungKut(:,1), rungKut(:,3), 'ro')

the resulting graph is displayed below.

Note that for data that is imprecisely aligned (i.e. not in the form of a matrix), the fscanf function can be used from within matlab. The fscanf function works in conjunction with the fopen and fclose functions and is used in a manner similar to their C counterparts. For more information, refer to the online help for these functions.
[Missing image]
Figure 7. Plot of Data Contained in a Data File

4.7. 3D Plots

MATLAB provides several functions for the visualization of 3-dimensional data. In Table 3 below, a list of these functions along witha brief description is presented.

MATLAB 3D Plot Functions

Function	Description
plot3	plots lines and points in 3D
contour, countour3	creates contour plots
pcolor	draws a rectangular array of cells whose colors are determined by matrix elements
image	draws a matrix as an image by mapping the elements of the matrix to the current color map
mesh, meshc, meshm	creates 3D perspective plots of matrix elements displayed as heights above an underlying plane.
surf, surfc, surfl	a combination of the mesh and contour functions
fill3	creates a 3D polygon and fills it with solid or interpolated colors.

Lines

The plot3 function works in a similar manner to the plot function. It requires an additional parameter that specifies the Z-coordinates of the points to be plotted. For example, the statements:

>>t = 0:pi/50:10*pi;
>>plot3(1.0./(1.0 + tan(t)), sin(t), t)

produces the following plot:
[Missing image]
Figure 8. Arbitrary 3D Plot

If the arguments to plot3d are matrices of the same size, then the columns of the 3 matrices are used as input data for the line plots.

4.7.2. Surfaces

The mesh function is used to generate surface plots of data. This function is particularly useful in visualizing large matrices. When mesh is invoked with a single matrix arguent, the values contained in the matrix are used as Z-axis data and the indices of the matrix are used as X- and Y-axes data.X- and Y-axes data can be explicitly specified by passing additional parameters to the mesh function. Note that if X- and Y-axes data are to be explicitly specified, then the matrices containing the X and Y data have to preceed the Z data matrix in the call to mesh. The X- and Y-axes data may also be explicitly passed to the mesh function using vectors instead of matrices. In such a case, the vectors are repeated to form a matrix of the correct size - in other words the Z data is plotted over a regularly spaced rectanguler grid. The order in which the vectors are utilized in the plot is:

(x(j), y(i), z(i,j))

The following examples illustrate these concepts.

>> x=-10:.5:10;
>> y=x;
>> [X Y]=meshgrid(x,y);
>> R=sqrt(X.^2+Y.^2) + eps;
>> Z=sin(R)./R;
>> mesh(Z);

The meshgrid function generates 2 matrices from the x and y vectors. The rows of the resulting matrix X are formed by repeating the vector x and the columns of the matrix Y are formed by repeating the vector y. The result of the mesh function is displayed below.
[Missing image]
Figure 9. Sombrero Plot

Note in the preceeding example that the X- and Y-axes range from 0 to 40 and not from -10 to 10 which is the range of the X and Y datasets. The following example illustrates the explicit specification of X- and Y-axes data:

>> x=-5:0.5:5;
>> y=x;
>> [X Y]=meshgrid(x,y);
>> R=sqrt(25.0 - X.^2 - Y.^2);
>> for i = 1:21           
      for j = 1:21           
         if imag(R(i,j)) ~= 0
            R(i,j) = 0+0i;   
         end
      end
   end
>> mesh(X,Y,R);

This example also illustrates the use of MATLAB programming constructs to selectively modify matrix elements. The do loops iterate over all elements of the matrix R and selectively replaces its imaginary elements with the value 0. The result of the call to mesh is displayed below.
[Missing image]
Figure 10. Hemisphere

4.8. Combining Multiple Plots on a Single Page

You can combine multiple plots on the same page by using calls to the subplot(a,b,x) command. This command breaks the figure window into a a-by-b matrix of small subplots and selects the x-th subplot for the current plot. The subplots are numbered starting at 1 and increasing along rows to the value mxn at bottom right of the matrix.

>> [X Y Z]=cylinder(4*cos(t));
>> subplot(2,2,1);
>> mesh(X);
>> subplot(2,2,2);
>> mesh(Y)
>> x=-2*pi:0.1:2*pi;
>> subplot(2,2,3);
>> plot(x,sin(x));
>> subplot(2,2,4);
>> plot(x, log(x);

The result of the sequence of commands above is displayed below.

[Missing image]
Figure 11. Multiple Plots on a Single Chart

4.1. - 2D Plots

4.6. - Plotting Data Contained in a Data File

4.7. - 3D Plots

4.7.1. - Lines

4.7.2. - Surfaces

4.8. - Combining Multiple Plots on a Single Page

Linear Discriminant Analysis Introduction

Introduction

Linear Discriminant Analysis (LDA) is a method to discriminate between two or more groups of samples. In order to develop a classifier based on LDA, you have to perform the following steps:
 

definition of groups

definition of discriminating function

estimation of discriminating function

test of discriminating function

application

Definition of groups:
The groups to be discriminated can be defined either naturally by the problem under investigation, or by some preceding analysis, such as a cluster analysis. The number of groups is not restricted to two, although the discrimination between two groups is the most common approach. Note that the number of groups must not exceed the number of variables describing the data set. Another prerequisite is that the groups have the same covariance structure (i.e. they must be comparable).
 
Definition of discriminating function:
In principle, any mathematical function may be used as a discriminating function. In case of the LDA, a linear function of the form

y = a₀ + a₁x₁ + a₂x₂ + ..... + a_nx_n is used, with x_i being the variables describing the data set. The parameters a_ihave to be determined in such a way that the discrimination between the groups is best. Note that this linear discriminating function is formally equivalent to the multiple linear regression. In fact, one can directly use MLR if the response variable y is replaced by the weighted class numbers c₁ and c₂:

c₁ = n₂/(n₁+n₂)    and    c₂ = - n₁/(n₁+n₂) In order to get a better understanding of the working of  the discriminating function, start the following .
 

Estimation of the parameters of the discriminating function: As you have seen in the interactive example above, there is only one direction of the discriminating line which yields the best separation results. The determination of the coefficients of the discriminating function is quite simple. In principle, the discriminating function is formed in such a way that the separation (=distance) between the groups is maximized, and the distance within the groups is minimized.
 
Test of the discriminating function
When the discriminating function is parametrized, it has to be tested either by using an independent set of test data, or by performing cross-validation. In both cases, the results of the test set should be comparable to the training data.
 
Application
Discriminant analysis can be used to perform either analysis or classification:
 

• Analysis: How can the material be interpreted? Which variables contribute most to the difference?
• Classification: Given that a discriminating function can be found which provides satisfactory separation, this function can be used to classify unknown objects.

source: http://www.vias.org/tmdatanaleng/cc_lda_intro.html

QRS Complex Detection and ECG Signal Processing

% QRS Detection Example
% shows the effect of each filter according to Pan-Tompkins algorithm.
% Note that, the decision  algorithm is different then the mentioned algorithm.

% by Faruk UYSAL

clear all
close all

x1 = load('ecg3.dat'); % load the ECG signal from the file
fs = 200;              % Sampling rate
N = length (x1);       % Signal length
t = [0:N-1]/fs;        % time index


figure(1)
subplot(2,1,1)
plot(t,x1)
xlabel('second');ylabel('Volts');title('Input ECG Signal')
subplot(2,1,2)
plot(t(200:600),x1(200:600))
xlabel('second');ylabel('Volts');title('Input ECG Signal 1-3 second')
xlim([1 3])

Cancellation DC drift and normalization

x1 = x1 - mean (x1 );    % cancel DC conponents
x1 = x1/ max( abs(x1 )); % normalize to one

figure(2)
subplot(2,1,1)
plot(t,x1)
xlabel('second');ylabel('Volts');title(' ECG Signal after cancellation DC drift and normalization')
subplot(2,1,2)
plot(t(200:600),x1(200:600))
xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second')
xlim([1 3])

Low Pass Filtering

% LPF (1-z^-6)^2/(1-z^-1)^2
b=[1 0 0 0 0 0 -2 0 0 0 0 0 1];
a=[1 -2 1];


h_LP=filter(b,a,[1 zeros(1,12)]); % transfer function of LPF

x2 = conv (x1 ,h_LP);
%x2 = x2 (6+[1: N]); %cancle delay
x2 = x2/ max( abs(x2 )); % normalize , for convenience .

figure(3)
subplot(2,1,1)
plot([0:length(x2)-1]/fs,x2)
xlabel('second');ylabel('Volts');title(' ECG Signal after LPF')
xlim([0 max(t)])
subplot(2,1,2)
plot(t(200:600),x2(200:600))
xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second')
xlim([1 3])

High Pass Filtering

% HPF = Allpass-(Lowpass) = z^-16-[(1-z^-32)/(1-z^-1)]
b = [-1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 32 -32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1];
a = [1 -1];

h_HP=filter(b,a,[1 zeros(1,32)]); % impulse response iof HPF

x3 = conv (x2 ,h_HP);
%x3 = x3 (16+[1: N]); %cancle delay
x3 = x3/ max( abs(x3 ));

figure(4)
subplot(2,1,1)
plot([0:length(x3)-1]/fs,x3)
xlabel('second');ylabel('Volts');title(' ECG Signal after HPF')
xlim([0 max(t)])
subplot(2,1,2)
plot(t(200:600),x3(200:600))
xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second')
xlim([1 3])

Derivative Filter

% Make impulse response
h = [-1 -2 0 2 1]/8;
% Apply filter
x4 = conv (x3 ,h);
x4 = x4 (2+[1: N]);
x4 = x4/ max( abs(x4 ));

figure(5)
subplot(2,1,1)
plot([0:length(x4)-1]/fs,x4)
xlabel('second');ylabel('Volts');title(' ECG Signal after Derivative')
subplot(2,1,2)
plot(t(200:600),x4(200:600))
xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second')
xlim([1 3])

Squaring

x5 = x4 .^2;
x5 = x5/ max( abs(x5 ));
figure(6)
subplot(2,1,1)
plot([0:length(x5)-1]/fs,x5)
xlabel('second');ylabel('Volts');title(' ECG Signal Squarting')
subplot(2,1,2)
plot(t(200:600),x5(200:600))
xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second')
xlim([1 3])

Moving Window Integration

% Make impulse response
h = ones (1 ,31)/31;
Delay = 15; % Delay in samples

% Apply filter
x6 = conv (x5 ,h);
x6 = x6 (15+[1: N]);
x6 = x6/ max( abs(x6 ));

figure(7)
subplot(2,1,1)
plot([0:length(x6)-1]/fs,x6)
xlabel('second');ylabel('Volts');title(' ECG Signal after Averaging')
subplot(2,1,2)
plot(t(200:600),x6(200:600))
xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second')
xlim([1 3])

Find QRS Points Which it is different than Pan-Tompkins algorithm

figure(7)
subplot(2,1,1)
max_h = max(x6);
thresh = mean (x6 );
poss_reg =(x6>thresh*max_h)';

figure (8)
subplot(2,1,1)
hold on
plot (t(200:600),x1(200:600)/max(x1))
box on
xlabel('second');ylabel('Integrated')
xlim([1 3])

subplot(2,1,2)
plot (t(200:600),x6(200:600)/max(x6))
xlabel('second');ylabel('Integrated')
xlim([1 3])

left = find(diff([0 poss_reg])==1);
right = find(diff([poss_reg 0])==-1);

left=left-(6+16);  % cancle delay because of LP and HP
right=right-(6+16);% cancle delay because of LP and HP

for i=1:length(left)
    [R_value(i) R_loc(i)] = max( x1(left(i):right(i)) );
    R_loc(i) = R_loc(i)-1+left(i); % add offset

    [Q_value(i) Q_loc(i)] = min( x1(left(i):R_loc(i)) );
    Q_loc(i) = Q_loc(i)-1+left(i); % add offset

    [S_value(i) S_loc(i)] = min( x1(left(i):right(i)) );
    S_loc(i) = S_loc(i)-1+left(i); % add offset

end

% there is no selective wave
Q_loc=Q_loc(find(Q_loc~=0));
R_loc=R_loc(find(R_loc~=0));
S_loc=S_loc(find(S_loc~=0));

figure
subplot(2,1,1)
title('ECG Signal with R points');
plot (t,x1/max(x1) , t(R_loc) ,R_value , 'r^', t(S_loc) ,S_value, '*',t(Q_loc) , Q_value, 'o');
legend('ECG','R','S','Q');
subplot(2,1,2)
plot (t,x1/max(x1) , t(R_loc) ,R_value , 'r^', t(S_loc) ,S_value, '*',t(Q_loc) , Q_value, 'o');
xlim([1 3])

Published with MATLAB® 7.10

source: http://farukuysal.net/QRSDetection.aspx

Using FFT to Obtain Simple Spectral Analysis Plots

Question:

How can you correctly scale the output of the FFT function to obtain a meaningful power versus frequency plot?

Answer:

Assume x is a vector containing your data. A sample vector used with this technical note is a 200 Hz sinusoid signal.

% Sampling frequency 
Fs = 1024;  


% Time vector of 1 second  
t = 0:1/Fs:1; 


% Create a sine wave of 200 Hz.
x =sin(2*pi*t*200);

First, you need to call the FFT function. For the fastest possible ffts, you will want to pad your data with enough zeros to make its length a power of 2. The built-in FFT function does this for you automatically, if you give a second argument specifying the overall length of the fft, as demonstrated below:

% Use next highest power of 2 greater than or equal to length(x) to calculate fft

nfft = 2^(nextpow2(length(x)));


% Take fft, padding with zeros so that length(fftx) is equal to nfft 

fftx = fft(x,nfft);

If nfft is even (which it will be, if you use the above two commands above), then the magnitude of the fft will be symmetric, such that the first (1+nfft/2) points are unique, and the rest are symmetrically redundant. The DC component of x is fftx(1) , and fftx(1+nfft/2)> is the Nyquist frequency component of x. If nfft is odd, however, the Nyquist frequency component is not evaluated, and the number of unique points is (nfft+1)/2 . This can be generalized for both cases to ceil((nfft+1)/2) .

% Calculate the number of unique points

NumUniquePts = ceil((nfft+1)/2);

% FFT is symmetric, throw away second half

fftx = fftx(1:NumUniquePts);

Next, calculate the magnitude of the fft:

% Take the magnitude of fft of x

mx = abs(fftx);

Consider the fact that MATLAB does not scale the output of fft by the length of the input:

% Scale the fft so that it is not a function of the length of x

mx = mx/length(x);


% Now, take the square of the magnitude of fft of x which has been scaled properly.
mx = mx.^2; 


% Since we dropped half the FFT, we multiply mx by 2 to keep the same energy.
% The DC component and Nyquist component, if it exists, are unique and should not be multiplied by 2. 

if rem(nfft, 2) % odd nfft excludes Nyquist point 
  mx(2:end) = mx(2:end)*2;
else
  mx(2:end -1) = mx(2:end -1)*2;
end

Now, create a frequency vector:

% This is an evenly spaced frequency vector with NumUniquePts points.


f = (0:NumUniquePts-1)*Fs/nfft;

Finally, generate the plot with a title and axis labels.

% Generate the plot, title and labels.
plot(f,mx);
title('Power Spectrum of a 200Hz Sine Wave');
xlabel('Frequency (Hz)'); 
ylabel('Power');

Bringing this all together, you get the following MATLAB file:

% Sampling frequency 

Fs = 1024; 


% Time vector of 1 second 

t = 0:1/Fs:1; 


% Create a sine wave of 200 Hz.

x = sin(2*pi*t*200); 


% Use next highest power of 2 greater than or equal to length(x) to calculate FFT.
nfft= 2^(nextpow2(length(x))); 


% Take fft, padding with zeros so that length(fftx) is equal to nfft 

fftx = fft(x,nfft); 


% Calculate the numberof unique points
NumUniquePts = ceil((nfft+1)/2); 


% FFT is symmetric, throw away second half 

fftx = fftx(1:NumUniquePts); 


% Take the magnitude of fft of x and scale the fft so that it is not a function of the length of x

mx = abs(fftx)/length(x); 


% Take the square of the magnitude of fft of x. 

mx = mx.^2; 


% Since we dropped half the FFT, we multiply mx by 2 to keep the same energy.
% The DC component and Nyquist component, if it exists, are unique and should not be multiplied by 2.


if rem(nfft, 2) % odd nfft excludes Nyquist point
  mx(2:end) = mx(2:end)*2;
else
  mx(2:end -1) = mx(2:end -1)*2;
end


% This is an evenly spaced frequency vector with NumUniquePts points. 

f = (0:NumUniquePts-1)*Fs/nfft; 


% Generate the plot, title and labels. 

plot(f,mx); 
title('Power Spectrum of a 200Hz Sine Wave'); 
xlabel('Frequency (Hz)'); 
ylabel('Power'); 

The resulting plot looks like the following:

The Signal Processing Toolbox 6.2 adds a new spectrum object. Spectrum objects contain parameter information for a particular spectral estimation method (e.g., spectrum.welch). This object provides an improved way to view and manipulate spectral estimation parameters. See the spectrum reference page (Spectrum Objects) and the associated estimation method reference pages (Spectral Estimation Method) for more information. The object has methods to evaluate the power spectral density for parametric and non-parametric (conventional) techniques and the mean-square spectrum for the non-parametric techniques. For subspace spectral estimation techniques (MUSIC and EIG), however, the object has methods to compute the pseudospectrum. As an example for using the spectrum object, refer to the following demo: Measuring the Power of Deterministic Periodic Signals.

source: http://www.mathworks.com/support/tech-notes/1700/1702.html

FFT and Spectral Leakage

at Thursday, January 20, 2011 Posted by mathuranathan

Introduction:

Frequency Transform is  used to study a signal's frequency domain characteristics. When using FFT to study the frequency domain characteristics of a signal, there are two limits : 1) The detectability of a small signal in the presence of a larger one ; 2) frequency resolution - which distinguishes two different frequencies.

In practice, the measured signals are limited in time and the FFT calculates the frequency transform over a certain number of discreet frequencies called bins.

Spectral Leakage:

In reality, signals are of time-limited nature and nothing can be known about the signal beyond the measured interval. For example, if the measurement of a never ending continuous train of sinusoidal wave is of interest, at some point of time we need to terminate our measurement to do further analysis. The limit on the time is also posed by limitations of the measurement system itself (like buffer size) besides other factors. Some assumptions have to made about the nature of the signal outside the measured interval. Fourier Transforms implicitly assumes that the signal essentially repeats itself after the measured interval.

Figure 1 illustrates the scenario in which a continuous train of sinusoidal signal is observed over a finite interval of time ("measured signal"). As discussed, the FFT assumes the signal to be continuous (conceptually, it does this by juxtaposing the measured signal repetitively). Observe the glitches in the assumed signal. These glitches are the manifestations of the measurement time relative to the frequency of the actual signal. If measurement time is an integral multiple of the rate of the actual signal (i.e. the inverse of the frequency of the signal), then no glitch will be observed in the assumed signal. In Figure 1, the measurement time is purposefully made to be a non-integral multiple of the actual signal rate. These sharp discontinuities will spread out in the frequency domain. This is called spectral leakage.

Figure 1: Impact of observation interval on FFT

Lets visualize the concept of spectral leakage by taking a pure sinusoidal signal as an example. Here we consider a sinusoid of 7Hz frequency (7 cycles in 1 second) and sampling it with a sampling frequency Fs=100Hz. We observe the signal for 100 seconds (700 cycles in total) and take the FFT of the observed signal. Figure 2 illustrates the frequency spectrum of the observed/measured signal. Essentially, the frequency spectrum contains a distinct peak at 7Hz.

Figure 2: Frequency Spectrum of a 7 Hz sinusoid observed for 100 seconds

Next, the measurement window is shrunk to 1 second. Now the FFT is taken on this observed signal. Figure 3 illustrates the frequency spectrum of the observed signal (which has 7 cycles only). The frequency spectrum contains some spectral leakage because of limited observation interval.

 Figure 3 : A 7 Hz sinusoid observed for 1 second and its Frequency Spectrum

In the previous illustrations (Figure 2 and 3), the observation time interval contained an integral number of sinewave cycles (i.e. in Figure 2 the observation interval contained exactly 700 cycles of sinusoid and in Figure 3 the observation interval contained exactly 7 cycles of sinusoid in the time domain).

For the next illustration the measurement time interval is adjusted in such a way that number of cycles in the observation window is no longer an integer. In Figure 4, the signal is observed for 1.4 seconds, which implies that there are 9.8 cycles. Now the observed signal does not end at zero amplitude at t=1.4 seconds. This scenario gives rise to gliches in the FFT's assumed signal (which constructs a periodic signal from the observed signal) and obviously results in more spectral leakage (Compare the Frequency spectrum in Figure 3 and Figure 4).
 

Figure 4 : A 7 Hz sinusoid observed for 1.4 seconds and its Frequency Spectrum.

 The effect of spectral leakage may be lessened if the observed signal does not contain any discontinuity at the end of the measurement time (this scenario rarely occurs in any real application). Another scenario in which the spectral leakage can be reduced is by having a signal that gradually reduces to zero at the ends of the measurement time.( All the windows (like Hamming, Hanning, Bartlett, etc.., essentially attempts to do this). Such signal would have no discontinuity when made periodic and so does not suffer spectral leakages.

Conclusion:

Spectral leakage is not due to FFT but due to the finite observation time. Spectral Leakage gives rise to two problems : 1) The spectral component of the desired signal no longer contains the complete energy, rather it also contains the energy of adjacent components and noise, thereby reducing the Signal to noise ratio; 2) Spectral leakage from a larger signal component may significantly overshadow other smaller signals making them difficult to identify or detect. 

Windowing Techniques are used to mitigate the effects of spectral leakage and therefore the restriction of having a limited observation interval.

Matlab Code:

Here is a matlab code to simulate the charts given in this discussion

%Illustration of Spectral Leakage and its effects on Frequency Spectrum

%(FFT)

%Author: Mathuranathan Viswanathan

%http:gaussianwaves.blosgpot.com

%License - Creative Commons - cc-by-nc-sa

observationTime=1; %Input : Observation time interval change it to 100 and 1.4

%and see the effect on FFT

Fx=7; %Frequency of the sinusoid

Fs=100; %Sampling Frequency
t=0:1/Fs:observationTime;
x=1*sin(2*pi*Fx*t);
plot(t,x)
title('SineWave - Frequency = 7Hz');
xlabel('Times (s)');ylabel('Amplitude');

%Perform FFT

NFFX=2.^(ceil(log(length(x))/log(2)));
FFTX=fft(x,NFFX);%pad with zeros
NumUniquePts=ceil((NFFX+1)/2);
FFTX=FFTX(1:NumUniquePts);
MY=abs(FFTX);
MY=MY*2;
MY(1)=MY(1)/2;
MY(length(MY))=MY(length(MY))/2;
MY=MY/length(x);
f1=(0:NumUniquePts-1)*Fs/NFFX;

figure;
stem(f1,MY,'k');
title('FFT of the sine wave');
axis([0 50 0 1]);
xlabel('Frequency');ylabel('Amplitude');

Reference:

[1] Fredric J Harris, "On the use of windows for Harmonic Analysis with the Discrete Wavelet Transform", Proceedings of IEEE, Vol 66,No 1, January 1978 - click here

Recommended Readings:

 

Read more: http://gaussianwaves.blogspot.com/2011/01/fft-and-spectral-leakage.html#ixzz1ICb1Capm
Under Creative Commons License: Attribution Non-Commercial Share Alike
source:http://gaussianwaves.blogspot.com/2011/01/fft-and-spectral-leakage.html

MATLAB - Introductory FFT Tutorial

Introduction

In this tutorial, we will discuss how to use the fft (Fast Fourier Transform) command within MATLAB. The fft command is in itself pretty simple, but takes a little bit of getting used to in order to be used effectively.
When we represent a signal within matlab, we usually use two vectors, one for the x data, and one for the y data. The fft command only operates on the y-data (converting the y-data from the time domain into the frequency domain), so it’s up to the user to determine what the x-data in the frequency domain will be! This tutorial will show you how to define your x-axis so that your fft results are meaningful. In addition, it will show you how to obtain a two-sided spectrum as well as a positive frequency spectrum for a given signal.

A Simple Example

1. Let’s start off with a simple cosine wave, written in the following manner:
   
   
2. Next, let’s generate this curve within matlab using the following commands:
   
   

   fo = 4;   %frequency of the sine wave
   Fs = 100; %sampling rate
   Ts = 1/Fs; %sampling time interval
   t = 0:Ts:1-Ts; %sampling period
   n = length(t); %number of samples
   y = 2*sin(2*pi*fo*t); %the sine curve
   
   %plot the cosine curve in the time domain
   sinePlot = figure;
   plot(t,y)
   xlabel('time (seconds)')
   ylabel('y(t)')
   title('Sample Sine Wave')
   grid

   Here’s what we get:
   
   
3. When we take the fft of this curve, we would ideally expect to get the following spectrum in the frequency domain (based on fourier theory, we expect to see one peak of amplitude 1 at -4 Hz, and another peak of amplitude 1 at +4 Hz):
   
   
4. There is also a phase component, but we’ll discuss that in a future tutorial.
   

Using Matlab’s FFT Command

So now that we know what to expect, let’s use MATLAB’s built in fft command to try to recreate the frequency spectrum:

%plot the frequency spectrum using the MATLAB fft command
matlabFFT = figure;  %create a new figure
YfreqDomain = fft(y); %take the fft of our sin wave, y(t)

stem(abs(YfreqDomain));  %use abs command to get the magnitude
%similary, we would use angle command to get the phase plot!
%we'll discuss phase in another post though!

xlabel('Sample Number')
ylabel('Amplitude')
title('Using the Matlab fft command')
grid
axis([0,100,0,120])

This doesn’t quite look like what we predicted above. If you notice, there are a couple of things that are missing.

• The x-axis gives us no information on the frequency. How can we tell that the peaks are in the right place?
• The amplitude is all the way up to 100
• The spectrum is not centered around zero

A Custom Function for fft to Obtain a Two-Sided Spectrum

Here is a helpful function I learned at Harvey Mudd College which will simplify the process of plotting a two-sided spectrum. Copy this code into an m-file and save it.

function [X,freq]=centeredFFT(x,Fs)
%this is a custom function that helps in plotting the two-sided spectrum
%x is the signal that is to be transformed
%Fs is the sampling rate

N=length(x);

%this part of the code generates that frequency axis
if mod(N,2)==0
    k=-N/2:N/2-1; % N even
else
    k=-(N-1)/2:(N-1)/2; % N odd
end
T=N/Fs;
freq=k/T;  %the frequency axis

%takes the fft of the signal, and adjusts the amplitude accordingly
X=fft(x)/N; % normalize the data
X=fftshift(X); %shifts the fft data so that it is centered

This is a relatively simple function to use. The function outputs the correct frequency range and the transformed signal. It takes in as input the signal to be transformed, and the sampling rate.
Let’s use the sine wave from above and do a quick example (Remember to set the Matlab directory to the location where you saved the previous m-file). Now, copy and paste these commands into the Matlab command prompt.

[YfreqDomain,frequencyRange] = centeredFFT(y,Fs);
centeredFFT = figure;

%remember to take the abs of YfreqDomain to get the magnitude!
stem(frequencyRange,abs(YfreqDomain));
xlabel('Freq (Hz)')
ylabel('Amplitude')
title('Using the centeredFFT function')
grid
axis([-6,6,0,1.5])

Here’s what you should see:

As you can see, this plot is basically identical to what we would expect! We get peaks at both -4 Hz and +4 Hz, and the amplitude of the peaks are 1.

Redundant Information in the FFT

As you can see from the plots above, the information within the frequency spectrum is entirely symmetric. Thus, we only need one side of the spectrum. In general, the positive side of the spectrum is used, while the negative side is ignored. So let’s adjust out function above so that we only get the positive frequencies.

A Custom Function for fft to Obtain only the Positive Frequencies

The following function is a modification of the above function, and will help you plot only the positive frequencies of the spectrum.

function [X,freq]=positiveFFT(x,Fs)
N=length(x); %get the number of points
k=0:N-1;     %create a vector from 0 to N-1
T=N/Fs;      %get the frequency interval
freq=k/T;    %create the frequency range
X=fft(x)/N; % normalize the data

%only want the first half of the FFT, since it is redundant
cutOff = ceil(N/2);

%take only the first half of the spectrum
X = X(1:cutOff);
freq = freq(1:cutOff);

Once again, let’s use the same sine wave and put it through this function. Copy and paste the following code into the Matlab command prompt.

[YfreqDomain,frequencyRange] = positiveFFT(y,Fs);
positiveFFT = figure;
stem(frequencyRange,abs(YfreqDomain));
set(positiveFFT,'Position',[500,500,500,300])
xlabel('Freq (Hz)')
ylabel('Amplitude')
title('Using the positiveFFT function')
grid
axis([0,20,0,1.5])

Here’s what you should get:

These two functions are very useful, and I still use them all the time!

Power of 2

The fft command within Matlab allows you to specify how many data points are in the transform. The Matlab documentation recommends that a power of 2 be used for optimal computation time. In my experience, there really isn’t a need to specify N as a power of 2. By using the next greatest power of 2, the fft command pads the original signal with zeros and proceeds to do a FFT on the signal. I’ve done some quick runs using fft with N as a power of 2 and N not as a power of 2, and the speed difference was neglible. In some cases, a 120 point FFT took LESS time than a 128 point FFT in some of my runs.
I don’t know exactly how the Matlab fft works, but I believe that it internally pads the signal with zeroes to the next greatest power of 2, performs the fft , then spits out an answer without the padded zeros.
I highly encourage anyone with greater knowledge to shed some light on this topic!

Recap and Future Topics

In this post, we talked primarily about how to use the fft command to create a frequency spectrum from a given signal. Two important things we did were to appropriately define the frequency axis, adjust the amplitude, and to view the spectrum as a two-sided, or one-sided spectrum.
In this post, we used a very simple example signal that was very well behaved. Unfortunately, this will rarely be the case. In the next upcoming posts, we will discuss zero-padding, windowing, filtering, and other techniques that will help you interpret and analyze the frequency spectrum of various signals.
Another thing that was not covered in this post, but is of imperative importance, is the phase. Stay tuned, and be ready for more material in the future!

Download source files and Further Reading

You can download the source files here.
Brush up your Fourier by reading about the theory and background at these links:
http://www.complextoreal.com/chapters/fft1.pdf
http://www.dspguide.com/ch8.htm
This is the end of the tutorial.

Reference:  http://blinkdagger.com/matlab/matlab-introductory-fft-tutorial/