Contents

% /home/dimarzio/Documents/working/12229/lectures/ftransform.m
% Sat Nov  4 14:39:30 2017
% Chuck DiMarzio, Northeastern University
%

Fourier transform

t=[0:4096-1]*10e-3/4096;  % Time axis
T=1024*10e-3/4096; % Period of a periodic function
f0=1/T
omega0=2*pi*f0;
% Sampling Frequency
dt=t(2)-t(1)
fsample=1/dt

f0 =
   400
dt =
   2.4414e-06
fsample =
      409600

Notation:

lower case v is time-domain voltage upper case v is frequency-domain voltage

Start with a sinusoidal wave

vm=5;phi=30*pi/180;v=vm*cos(omega0*t+phi);
figure;plot(t,v,'b-');grid on;
title('vm=5;phi=30*pi/180;v=vm*cos(omega0*t+phi)')
xlabel('t, Time, sec');ylabel('v, Volts');
%  Fourier transform
V=fft(v);
figure;cxplot(1:length(t),V);grid on;
title('V=fft(v)');
xlabel('?');ylabel('?');
% Plot the first and last 10 points
figure;cxplot(1:10,V(1:10));grid on;
title('V=fft(v) Zoom low end');
xlabel('?');ylabel('?');
figure;cxplot(4087:4096,V(4087:4096));
grid on;
title('V=fft(v) Zoom high end');
xlabel('?');ylabel('?');

Use a frequency axis and scale the voltage

V=V/4096*max(t);
f=fftaxis(t); % needs fftaxis.m by Chuck DiMarzio
% Plot the first and last 10 points
figure;cxplot(f(1:10),V(1:10));grid on;
xlabel('f, Frequency, Hz');ylabel('V, Volts/Hz');
title('V=fft(v) zoom low end with f axis');
figure;cxplot(f(4087:4096)/1000,V(4087:4096));
grid on;
title('V=fft(v) zoom high end with f axis');
xlabel('F, Frequency, kHz');ylabel('V, Volts/Hz');

Deal with the periodicity of the FFT

f=fftaxisshift(f); % needs fftaxisshift.m by Chuck DiMarzio
V=fftshift(V); % Existing function in matlab
figure;cxplot(f,V);grid on;
xlabel('F, Frequency, kHz');ylabel('V, Volts/Hz');
title('V=fft(v) with fftshift');
%Plot the interesting region
test=find(abs(f)<2*f0);
figure;cxplot(f(test),V(test));grid on;
title('V=fft(v) with fftshift zoom to center');
xlabel('F, Frequency, kHz');ylabel('V, Volts/Hz');
disp('What is the amplitude?');
disp('There are only 2 points, so inverse FT is quite easy.');
disp('Multiply by df and exp(1j*2*pi*f*t) where f = +/- f0.');
Vmax=max(abs(V))

What is the amplitude?
There are only 2 points, so inverse FT is quite easy.
Multiply by df and exp(1j*2*pi*f*t) where f = +/- f0.
Vmax =
    0.0250

Now let's do a square wave with DC bias

vm=5;phi=0*pi/180;v=double(vm*cos(omega0*t+phi)>0);
figure;plot(t,v,'b-');grid on;
title('vm=5;phi=0*pi/180;v=double(vm*cos(omega0*t+phi)>0)');
xlabel('t, Time, sec');ylabel('v, Volts');
%  Fourier transform
V=fftshift(fft(v)/4096*max(t));
test=find(abs(f)<20*f0);
figure;cxplot(f(test),V(test));grid on;
title('FT of Square Wave');
xlabel('F, Frequency, Hz');ylabel('V, Volts/Hz');

Now a pulse

v=zeros(size(t));
v(1:256)=1;
figure;plot(t,v,'b-');grid on;
title('Pulse');
xlabel('t, Time, sec');ylabel('v, Volts');
%  Fourier transform
V=fftshift(fft(v)/4096*max(t));
test=find(abs(f)<20*f0);
figure;cxplot(f(test),V(test));grid on;
title('FT of Pulse zoom to center');
xlabel('F, Frequency, Hz');ylabel('V, Volts/Hz');

OK. time for an interesting sinusoid

Don't try to understand this until you know about convolution

T1=900.354*10e-3/4096; % Not a power of 2
f1=1/T1
omega1=2*pi*f1;
vm1=5;phi1=0*pi/180;v=vm1*cos(omega1*t+phi1);
figure;plot(t,v,'b-');grid on;
title('Cosine Wave');
xlabel('t, Time, sec');ylabel('v, Volts');
%  Fourier transform
V=fftshift(fft(v)/4096*max(t));
test=find(abs(f)<5*f0);
figure;cxplot(f(test),V(test));grid on;
title('FT Cosine Wave zoom to center');
xlabel('F, Frequency, Hz');ylabel('V, Volts/Hz');

f1 =
  454.9322

Truncate the sine wave

vt=v;
vt(find(t>2.95e-3))=0;
figure;plot(t,v,'b--',t,vt,'g-');grid on;
title('Truncated Cosine Wave');
xlabel('t, Time, sec');ylabel('v, Volts');
%  Fourier transform
Vt=fftshift(fft(vt)/4096*max(t));
test=find(abs(f)<5*f0);
figure;cxplot(f(test),Vt(test));grid on;
title('FT Truncated Cosine Wave zoom to center');
xlabel('F, Frequency, Hz');ylabel('V, Volts/Hz');

% Compare the Fourier Transform of a "time window"
w=ones(size(t));
w(find(t>2.95e-3))=0;
figure;plot(t,w,'m-');grid on;
title('Window Function');
xlabel('t, Time, sec');ylabel('Window');
%  Fourier transform
W=fftshift(fft(w)/4096*max(t));
test=find(abs(f)<5*f0);
figure;cxplot(f(test),W(test));grid on;
title('FT Window Function');
xlabel('F, Frequency, Hz');ylabel('FT of Window');

Now. let's add some noise

vn=vt+randn(size(t))*8;
figure;
plot(t,vn,'Color',0.7*[1,1,1]);grid on;
hold on;plot(t,v,'b-',t,vt,'g--');hold off;
title('Truncated Cosine with noise');
xlabel('t, Time, sec');ylabel('v, Volts');
%  Fourier transform
Vn=fftshift(fft(vn)/4096*max(t));
test=find(abs(f)<5*f0);
figure;cxplot(f(test),Vn(test));grid on;
title('FT Truncated Cosine with noise');
xlabel('F, Frequency, Hz');ylabel('V, Volts/Hz');

figure;plot(f(test),20*log10(abs(V(test))),'g-',...
	    f(test),20*log10(abs(Vt(test))),'m-',...
	    f(test),20*log10(abs(Vn(test))),'k-');grid on;
title('Compare Ft Cosines zoom to center');
xlabel('F, Frequency, Hz');ylabel('V, Volts/Hz');
legend('Full Signal','Truncated','Truncated+Noise');

% Do it in dB over a wider band
test=find(abs(f)<50*f0);
figure;plot(f(test)/1000,20*log10(abs(V(test))),'g-',...
	    f(test)/1000,20*log10(abs(Vt(test))),'m-',...
	    f(test)/1000,20*log10(abs(Vn(test))),'k-');grid on;
title('Compare Ft Cosines wideband');
xlabel('F, Frequency, kHz');ylabel('V, Volts/Hz');
legend('Full Signal','Truncated','Truncated+Noise');

Published with MATLAB® R2017a