function bracewell % Illustrates some properties of the Fourier transform % Calculates the Fourier transform of a signal, and plots its real and % imaginary parts, and its squared absolute value divided by the sampling % frequency (the one-sided spectral density). % % Last modified by fjsimons-at-alum.mit.edu, 6.10.2005 % Length of the sequence, in physical units, say seconds tmin=0; tmax=10*pi; noil=4; % Number of time points nsamp=2^9; % Number of fft points: equal to nsamp nfft=nsamp; % Discrete wavelengths present in the signal wavl=[0.5 0.2]; % Amplitudes of these components ampl=[1 1]; % Total time ttim=tmax-tmin; % Sampling frequency Fs=nsamp/ttim; % Time-domain signal t=linspace(tmin,tmax,nsamp); % Construct the time-domain signal s=zeros(size(t)); % Individual sine harmonic components for index=1:length(wavl) s=s+ampl(index)*sin(t*2*pi/wavl(index)); end % Add noise s=s+randn(1,length(t))/noil; % Frequency-domain signal: the Fourier transform F=fft(s,nfft); % The first and the last elements are real, the others complex % Note that the zero-frequency equals the sum of the data disp(sprintf('mean(s)= %8.3f ; F(1)/N= %8.3f ; diff= %8.3e',... mean(s),F(1)/length(s),abs(mean(s)-F(1)/length(s)))) % And Parseval's relation that the variance of the data is the variance disp(sprintf('sqnorm(s)= %8.3f ; sqnorm(F)/N= %8.3f ; diff= %8.3e',... s*s',F*F'/length(s),abs(s*s'-F*F'/length(s)))) % Frequency axis [f,selekt,l]=fftaxis1D(s,nfft,ttim); % Select the nonredundant components of the "spectrum" P=abs(F(selekt)).^2; % Thus, nfft=nsamp, but get imaginary and real, only half of it is unique, % thus length(real(F(selekt)))+length(imag(F(selekt)))-2=nsamp again % But don't forget to double the formerly redundant ones P=[P(1) 2*P(2:end-1) P(end)]; % 'Cause if abs(sum(P)-F*F')>1e-8 error('You are in trouble') end % Spectral density S=P/Fs; % What is the integral of S over all these frequencies? ig=trapeze(f,S); disp(sprintf('sqnorm(s)= %8.3f ; int(S(f))= %8.3f ; diff= %5.3f',... s*s',ig,abs(s*s'-ig))) % The accuracy here depends on the integration routine, rather clf ah=krijetem(subnum(2,2)); % Plot these various things and annotate them axes(ah(1)) p=plot(t,real(s)); tl(1)=title('Time Domain, Real'); axis([0 tmax minmax(real(s))*1.1]) grid on xl(1)=xlabel('t, time (s)'); yl(1)=ylabel('signal amplitude (a)'); axes(ah(2)) plot(f,S) %tl(2)=title('Spectral density'); axis([minmax(f) minmax(S)*1.1]) %movev(tl(2),range(ylim)/15) grid on xl(2)=xlabel('f, frequency (s^{-1})'); yl(2)=ylabel('spectral density (a^2s)'); axes(ah(3)) plot(real(fftshift(F))) tl(3)=title('Fourier transform'); axis([1 length(F) minmax(real(F))*1.1]) grid on xl(3)=xlabel('frequency number'); yl(4)=ylabel('real part (a)'); xval=100; xels=unique([1 xval:xval:nfft nfft]); xels=[xels(diff(xels)>xval/2) xels(end)]; set(gca,'xtick',xels) axes(ah(4)) plot(imag(fftshift(F))) tl(4)=title('Fourier transform'); axis([1 length(F) minmax(imag(F))*1.1]) grid on yl(4)=ylabel('Imaginary part (a)'); xl(4)=xlabel('frequency number'); set(gca,'xtick',xels) fig2print(gcf,'portrait') % Wavelength axis on second plot xtra=xtraxis1dlin(ah(2),[1/8 1/6 1/4 1/2 1 floor(ttim)],... {num2str(floor(ttim)) '1' '1/2' '1/4' '1/6' '1/8'},1); axes(xtra) xlx=xlabel('\lambda, Wavelength (s)'); movev(xlx,-range(ylim)/200) % Looks bad but prints good longticks([ah xtra]) figdisp