% EE3220 - Tuesday 1/26/2010 - Dr. Durant

% MATLAB DFT/IDFT example with phase



% In this file we construct a 7-sample signal with 3 components, including one 

% that is phase shifted, have MATLAB calculate its DFT, speculate what its DFT 

% should be based on our knowledge of the 2-sided phase spectrum, and then verify 

% that the two are the same.



% Let's start by constructing a 7-sample signal that has 3 components

% * 0th harmonic (DC)          , magnitude -5

% * 1st harmonic (Fundamental) , magnitude 3, cosine with phase lag of 45 degrees

% * 2nd harmonic (1st overtone), nothing

% * 3rd harmonic (2nd overtone), magnitude 1, cosine, 



% Note that the 3nd harmonic will go through 3 complete cycles (by definition),

% but that each cycle does NOT line up with an integer number of samples.

% Each cycle will take N/3 = 7/3 = 2.333 samples.  This is perfectly fine.



format compact

N = 7;

n = 0 : (N-1);



% Here are the 3 components by frequency

xDC  = -5 * cos (2*pi * 0 * n / N)

x1st =  3 * cosd(360  * 1 * n / N - 45)

x3rd =  1 * cos (2*pi * 3 * n / N)



% Add the 3 components into a single signal

x = xDC + x1st + x3rd



% Plot with circled points, connecting lines, and a legend...

plot(n,xDC,'bo-', n,x1st,'go-', n,x3rd,'ro-', n,x,'ko-')

legend('DC','1st harmonic','3rd harmonic','sum')



c = fft(x);



% without going through the detailed math, we expect the following values

c0 = N * -5; % DC level

c1 = N/2 * 3 * exp(-j*pi/4) % delay = negative phase added

c2 = 0;

c3 = N/2 * 1

c4 = conj(c3); 

c5 = conj(c2);

c6 = conj(c1);

% Notes on c4+: c4 = conjugate(c[N-4]) = conjugate(c[3]); we're lowering the 

% frequency of the complex exponential by 2pi (no effect due to period), so we 

% arrive at the negative frequency, which we know has opposite phase due to 

% Euler's decomposition of the cosine function. So c3&c4 (c3 & c-3) form a pair 

% and the /2 divides the real signal energy between them.



% Now put them together in an array:

cc = [c0 c1 c2 c3 c4 c5 c6]



% If we interpreted everything above correctly, c, calculated using the 

% DFT summation formula by MATLAB, and cc, calculated using our understanding 

% of 2-sided spectra, should agree. Let's calculate the difference, 

% element-by-element:

diff = c - cc



% It looks like everything is pretty close to 0.  To get a single number to 

% evaluate how close it is, let's add up the squared magnitudes of the errors (so 

% we don't get out of phase errors appearing to cancel).

SSE_c = sum(abs(diff).^2)



% Mathematically, this is equivalent to taking the dot product of the complex error vector with itself...

SSE_c2 = dot(diff,diff)

% This is about 5 * 10^-29, which is virtually 0.  It is not exactly 0 due to rounding errors



% And, we can invert the transform to get back to the original data:

xx = ifft(cc)



% And, confirm that the reconstructed signal equals the original (error is ~ 0)...

SSE_x = sum((x-xx).^2)