aguidetothefft.m

FFT的MATLAB程序

%-------------------------------------------------------------------
% Daniele Disco
% 10.08.2004
%-------------------------------------------------------------------
% GENERAL DESCRIPTION
% In the following, temporal signals will be characterized by 
% lower letters and the Fourier Transform with capilats letters:
%           x(t) === FT ===> X(f).
%
% The analysis are related to  
%  o FFT normalization;
%  o Verify the symmetry properties of the FFT
%  o Use of the FFT to produce a delayed segnal:
%        x(t) === FFT ===> X(f);
%        X(f) exp[-j(2*pi*f*taun)] === IFFT ===> x(t - taun).
%        This algorithm is useful when the delay is not a multiple
%        of the sampling time.
%-------------------------------------------------------------------

%-------------------------------------------------------------------
%-------------------------------------------------------------------
%-------------------------------------------------------------------
% Symmetric Time Domain
%-------------------------------------------------------------------
%-------------------------------------------------------------------
%-------------------------------------------------------------------
% The signals:
%     x(t) = exp[- t^2 / ( 2 A^2 )].
% The analityc FT:
%     X(f) = T sqrt(2 pi) exp(- 2 pi^2 f^2 A^2)
% The bandwidth a 3dB
%     Bx = sqrt(log(2)) / (2 pi A );
%-------------------------------------------------------------------

%-------------------------------------------------------------------
% Reset
%-------------------------------------------------------------------
clear all
close all
do_print = 1;
clc
%-------------------------------------------------------------------

%-------------------------------------------------------------------
% Parameters definition
%-------------------------------------------------------------------
Bx = 10;
A = sqrt(log(2))/(2*pi*Bx);  
fs = 500;                      
dt = 1/fs; 
T=1;
t = -T/2:dt:T/2-dt;                  
df = 1/T;
Fmax = 1/2/dt;
f=-Fmax:df:Fmax-df;
taun = 0.055;
%-------------------------------------------------------------------

%-------------------------------------------------------------------
%  Numerical evaluations
%-------------------------------------------------------------------
x = exp(-t.^2/2/A^2);
Xan = A*sqrt(2*pi)*exp(-2*pi^2*f.^2*A^2);                 
X = dt * fftshift(fft(fftshift(x)));
Xbad = dt * fftshift(fft(x));
ifftx = fs * ifftshift(ifft(ifftshift(X)));
esp = exp(-j * 2 * pi * f * taun);
Xshift = X .*esp;
Xanshift = Xan .* esp;
ifftxshift = fs * ifftshift(ifft(ifftshift(Xshift)));
%-------------------------------------------------------------------


%-------------------------------------------------------------------
% Graphical output
%-------------------------------------------------------------------
% Fig.0
%-------------------------------------------------------------------
figure;
plot(f,real(Xbad))
xlabel('Frequency [Hz]');
ylabel('fft()');
axis([-50 50 -0.04 0.04])
grid;
if(do_print);print -depsc2 fig0.eps;end

%-------------------------------------------------------------------
% Fig.1
%-------------------------------------------------------------------
figure;
subplot(2,1,2);
plot(t,x,'k');
xlabel('Time [s]');
ylabel('x(t)');
axis([-0.1 0.1 0 1])
grid
subplot(2,1,1);
plot(f,Xan,'k');
xlabel('Frequency [Hz]');
ylabel('X(f)');
axis([-250 250 0 0.04])
grid;
if(do_print);print -depsc2 fig1.eps;end
%-------------------------------------------------------------------
% Fig.2
%-------------------------------------------------------------------
figure;
plot(f,20*log10(Xan/(A*sqrt(2*pi))),'k');hold on;
xlabel('Frequency [Hz]');
ylabel('|X(f)| [dB]');
plot([+Bx +Bx],[-10 -3],'r');
plot([-Bx -Bx],[-10 -3],'r');
plot([-Bx +Bx],[-3 -3],'r');
h=text(-1.7, -2.6, '2B_x','fontweight','bold','fontangle','italic');
axis([-50 50 -10 0]);
grid
if(do_print);print -depsc2 fig2.eps;end
%-------------------------------------------------------------------
% Fig.3
%-------------------------------------------------------------------
figure;
subplot(2,1,1)
plot(f,real(X),'k'); hold on
plot(f,Xan,'r--');
xlabel('Frequency [Hz]');
ylabel('X(f)');
legend('Num.','Ana.');
axis([-250 250 0 0.04])
grid
subplot(2,1,2)
plot(t,real(ifftx),'k'); hold on
plot(t,x,'r--');
xlabel('Time [s]');
ylabel('x(t)');
legend('Num.','Ana.');
axis([-0.1 0.1 0 1]);
grid
if(do_print);print -depsc2 fig3.eps;end
%-------------------------------------------------------------------
% Fig.4
%-------------------------------------------------------------------
figure;
subplot(2,1,1);
plot(f,real(Xshift),'k'); hold on
plot(f,real(Xanshift),'r--');
title('X(f) exp(-j2 \pi f \tau_{n})');
legend('Re[Num.]','Re[Ana.]');
grid;
hold off
subplot(2,1,2);
plot(f,imag(Xshift),'k'); hold on
plot(f,imag(Xanshift),'r--');
legend('Im[Num.]','Im[Ana.]');
grid;
xlabel('Frequency [Hz]');
if(do_print);print -depsc2 fig4.eps;end;
%-------------------------------------------------------------------
% Fig.5
%-------------------------------------------------------------------
figure;
plot(t(1),real(ifftxshift(1)),'r',t,x,'k--',t,real(ifftxshift),'k');hold on
xlabel('Time [s]');
ylabel('x(t - \tau_{n})');
plot([taun taun],[0 1],'r');
legend(['\tau_n = 0.055s = ' num2str(taun/dt,3) 'dt']);
grid;
axis([-0.1 0.3 0 1]);
if(do_print);print -depsc2 fig5.eps;end
%-------------------------------------------------------------------

%-------------------------------------------------------------------
%-------------------------------------------------------------------
%-------------------------------------------------------------------
% Asymmetric Time Domain
%-------------------------------------------------------------------
%-------------------------------------------------------------------
%-------------------------------------------------------------------
% The signals:
%     x(t) = A t exp(-A t)     t>=0
% The analityc FT:
%                    A
%     X(f) = ------------------
%            ( A + 2 j pi f)^2
% The bandwidth a 3dB
%     Bx = A sqrt(sqrt(2)-1) / (2 pi);
%-------------------------------------------------------------------

%-------------------------------------------------------------------
% Reset only memory
%-------------------------------------------------------------------
clear all
do_print = 1;
clc
ui = complex(0,1);
%-------------------------------------------------------------------

%-------------------------------------------------------------------
% Parameters definition
%-------------------------------------------------------------------
Bx = 5;
A = 2*pi*Bx/sqrt(sqrt(2)-1);
fs = 500;                      
dt = 1/fs;                      
T = 1;                       
t1 = 0:dt:T/2-dt;
t = -T/2:dt:T/2-dt;
L0 = length(t)-length(t1);
df = 1/T;
Fmax = 1/2/dt;
f=-Fmax:df:Fmax-df;
taun = 0.055;
%-------------------------------------------------------------------

%-------------------------------------------------------------------
%  Numerical evaluations
%-------------------------------------------------------------------
x = A*t1.*exp(-A*t1);
x = [zeros(1,L0) x];
Xan = A./(A+ui*2*pi.*f).^2;                  
X = dt * fftshift(fft(fftshift(x)));                     
ifftx = fs * ifftshift(ifft(ifftshift(X)));
esp = exp(-j * 2 * pi * f * taun);
Xshift = X .*esp;
Xanshift = Xan .* esp;
ifftxshift = fs * ifftshift(ifft(ifftshift(Xshift)));
%-------------------------------------------------------------------

%-------------------------------------------------------------------
% Graphical output
%-------------------------------------------------------------------
% Fig.6
%-------------------------------------------------------------------
figure;
subplot(2,1,2)
plot(t,x,'k');
xlabel('Time [s]');
ylabel('x(t)');
grid
subplot(2,1,1)
[AX,H1,H2] = plotyy(f,real(Xan),f,imag(Xan),'plot');
set(get(AX(1),'Ylabel'),'String','Re[X(f)]');
set(get(AX(2),'Ylabel'),'String','Im[X(f)]');
set(AX(1),'Ycolor','k');
set(AX(2),'Ycolor','r');
set(H1,'Color','k')
set(H2,'Color','r')
grid;
title('X(f)');
xlabel('Frequency [Hz]')
if(do_print);print -depsc2 fig6.eps;end
%-------------------------------------------------------------------
% Fig.7
%-------------------------------------------------------------------
figure;
plot(f,20*log10(abs(Xan)*A),'k');hold on;
xlabel('Frequency [Hz]');
ylabel('Amplitude [dB]');
plot([+Bx +Bx],[-10 -3],'r');
plot([-Bx -Bx],[-10 -3],'r');
plot([-Bx +Bx],[-3 -3],'r');
text(-1.7, -2.6, '2B_x','fontweight','bold','fontangle','italic');
axis([-50 50 -10 0]);
grid
if(do_print);print -depsc2 fig7.eps;end
%-------------------------------------------------------------------
% Fig.8
%-------------------------------------------------------------------
figure;
subplot(3,1,1);
plot(f,real(X),'k'); hold on
plot(f,real(Xan),'r--');
ylabel('Re[X(f)]');
legend('Num.','Ana.');
grid
subplot(3,1,2);
plot(f,imag(X),'k'); hold on
plot(f,imag(Xan),'r--');
xlabel('Frequency [Hz]');
ylabel('Im[X(f)]');
legend('Num.','Ana.');
grid
subplot(3,1,3);
plot(t,real(ifftx),'k'); hold on
plot(t,real(x),'r--');
xlabel('Time [s]');
ylabel('x(t)');
legend('Num.','Ana.');
grid
if(do_print);print -depsc2 fig8.eps;end
%-------------------------------------------------------------------
% Fig.9
%-------------------------------------------------------------------
figure;
subplot(2,1,1);
plot(f,real(Xshift),'k'); hold on
plot(f,real(Xanshift),'r--');
title('X(f) exp(-j2 \pi f \tau_{n})');
legend('Re[Num.]','Re[Ana.]');
grid;
hold off
subplot(2,1,2);
plot(f,imag(Xshift),'k'); hold on
plot(f,imag(Xanshift),'r--');
legend('Im[Num.]','Im[Ana.]');
grid;
xlabel('Frequency [Hz]');
if(do_print);print -depsc2 fig9.eps;end
%-------------------------------------------------------------------
% Fig.10
%-------------------------------------------------------------------
figure;
plot(t(1),real(ifftxshift(1)),'r',t,x,'k--',t,real(ifftxshift),'k');hold on
xlabel('Time [s]');
ylabel('x(t - \tau_{n})');
plot([taun taun],[0 1],'r');
legend(['\tau_n = 0.055s = ' num2str(taun/dt,3) 'dt']);
grid;
axis([-0.1 0.2 0 0.4]);
if(do_print);print -depsc2 fig10.eps;end
%-------------------------------------------------------------------

%-------------------------------------------------------------------
%-------------------------------------------------------------------
%-------------------------------------------------------------------
% Periodic function
%-------------------------------------------------------------------
%-------------------------------------------------------------------
%-------------------------------------------------------------------

%-------------------------------------------------------------------
% Reset only memory
%-------------------------------------------------------------------
clear all
clc
do_print = 1;
%-------------------------------------------------------------------

%-------------------------------------------------------------------
% Parameters definition
%-------------------------------------------------------------------
f0 = 50;
fs = 500;                      
dt = 1/fs; 
T=1;
t = -T/2:dt:T/2-dt;                  
df = 1/T;
Fmax = 1/2/dt;
f=-Fmax:df:Fmax-df;
%-------------------------------------------------------------------

%-------------------------------------------------------------------
%  Numerical evaluations
%-------------------------------------------------------------------
x = cos(2*pi*f0*t); 
X = fftshift(fft(fftshift(x)))/length(x);     
ifftx = fs * ifftshift(ifft(ifftshift(X)));
%-------------------------------------------------------------------

%-------------------------------------------------------------------
% Graphical output
%-------------------------------------------------------------------
% Fig.11
%-------------------------------------------------------------------
figure; 
plot(t,x,'k');
xlabel('Time [s]');
ylabel('x(t)');
axis([-0.1 0.1 -1.5 1.5])
grid
if(do_print);print -depsc2 fig11.eps;end
%-------------------------------------------------------------------
% Fig.12
%-------------------------------------------------------------------
figure;
plot(f,real(X),'r');
ylabel('X(f)');
xlabel('Frequency [Hz]');
axis([-150 150 0 0.6])
grid
if(do_print);print -depsc2 fig12.eps;end
%-------------------------------------------------------------------
% Fig.13
%-------------------------------------------------------------------
figure;
plot(t,real(ifftx))
grid
axis([-0.1 0.1 -1.5 1.5])
if(do_print);print -deps fig13.eps;end
ylabel('x(t)=FFT^{-1}[X(f)]');
xlabel('time [s]');
if(do_print);print -depsc2 fig13.eps;end