% Script 3 for Week 6
So we are done with Fourier series and on to the Discrete Fourier Transform.
First download and place the two main function M-files for DFT in the MATLAB working directory: DFTtest.m and RealDFT.m. The first one is really long, but many Bothans died to bring us this information. For convenience I restrict this second one to deal with an even number of data points. The first has no such restriction.
The Discrete Fourier transform (DFT) is nothing more than the conversion of ordinary coordinates in R^N to coordinates in a very particular new orthonormal basis of R^N or, sometimes, C^N.
In the situation of interest N will be large (sometimes VERY large) and consist of some kind of measured values over EQUAL TIME INCREMENTS. Think about pressure measurements taken every second for N seconds. A vector in R^N in this context is called a record and its components are called samples.
You may recall that if F is the coordinate vector that represents the record in the new basis and X is the original record and M is the matrix of transition, containing as its columns the new basis vectors in the usual basis, then X=M*F. That means F=M\X. That is all there is to it.
Unfortunately, the left division indicated takes O(N^3) flops to perform, and if N=44000 (this is not a randomly chosen value) then N^3 is about 10^(14), a nonstarter.
Using a related basis in C^N we can do way better. It seems like everything should be more complicated using complex numbers and complex vectors, but because of extreme symmetries inside this new M there is a way of heavily recycling parts of the calculation in M\X and this method is called Fast Fourier Transform, FFT.
So the FFT is a way of quickly calculating the DFT, which is what we really want.
It takes O(N*Log(N)) flops to calculate the FFT. In the N=44000 case that is about about 2 million flops and yields a speed-up of a factor greater than 100 million.
We can use DFT methods to estimate the frequency and amplitude of underlying sinusoidal components to data X through MATLAB directly (not by FFT) using my script by:
[df M] = RealDFT(X)
where df is the vector of coefficients and M is the matrix whose columns are the elemental periodic vectors.
We assume that the length of data vector X is even (for convenience) and we also assume that the data was taken at equal time increments (a requirement).
CHANGE SCALE IF NECESSARY SO THIS TIME INCREMENT IS 1. It is a one-step operation to go back, and doing it this way makes many formulas easier to look at. So X(1) is the data taken at time 0 and X(2) is the data taken at time 1 etc.
If length(X)=N the first N/2+1 entries of df (and columns of M) correspond
to cosine terms and the last N/2-1 correspond to sine terms.
Except for m=1 and m=N/2+1 the cosine term m will have the
same frequency as the sine term m+N/2
The cosine terms for m=1 and m=N/2+1 are special: the m=1 cosine term has period 0, so is constant. The N/2+1 cosine term has period 2 corresponding to the Nyquist frequency. The sine terms for both these would be identically 0, so they are not included.
y=[3;7;9;1;5;-2]
N=length(y)
[df M omegas W ff] = RealDFT(y);
df
M
% M*df=y so df=M\y.
df(1)*M(:,1)+df(2)*M(:,2)+df(5)*M(:,5)+df(3)*M(:,3)+df(6)*M(:,6)+df(4)*M(:,4)
% When you start with an even number of data points there will be
% two more cosine terms than sine terms.
% for m=2 to N/2 the vectors
df(m)*M(:,m)+df(m+N/2)*M(:,m+N/2)
df(2)*M(:,2)+df(2+N/2)*M(:,2+N/2)
corresponds to the combined sine and cosine terms with period N/(m-1) and frequency (m-1)/N.
The following is N times the same vectors but calculated using W with ff coefficients
ff(m)*W(:,m)+ff(N+2-m)*W(:,N+2-m)
These pairs start at front-back and work "inwards."
The squared norm
norm(df(m)*M(:,m)+df(m+N/2)*M(:,m+N/2)).^2
norm(df(2)*M(:,2)+df(2+N/2)*M(:,2+N/2)).^2
is called the power at that frequency and is useful to compare the relative importance of different frequencies in the signal. The first coefficient is not important for this purpose, corresponding to a constant background signal.
for m=2:(N/2)
p(m)=norm(df(m)*M(:,m)+df(m+N/2)*M(:,m+N/2)).^2;
end
p(1)=norm(df(1)*M(:,1)).^2;
p(N/2+1)=norm(df(N/2+1)*M(:,N/2+1)).^2
for m=2:(N/2)
q(m)=norm(ff(m)*W(:,m)+ff(N+2-m)*W(:,N+2-m)).^2;
end
q(1)=norm(ff(1)*W(:,1)).^2;
q(N/2+1)=norm(ff(N/2+1)*W(:,N/2+1)).^2
q=q/N^2
%*********************************%************************************
function z = elemental(n)
% This will plot the n+1 elemental periodic functions on [0,n].
% It illustrates one reason why aliasing occurs.
t=linspace(0,n,1000);
t1=0:n;
k=0;
figure(12)
clf
subplot(2,1,1)
hold on
while k