Matlab code for FFT book

Data files for Chapters 8-9

Root-raised-cosine shaped data

Nonlinear device output signal

Intuitive Guide to Fourier Analysis and Spectral Estimation
Matlab code for key figures and examples is given here.

If you find any issues with these programs, please email me.
Last edited 11/24/2017

Matlab files from the book are included here.

% Chapter 1

%% Figure 1.8 - Summation of sines and cosine waves

clf
time = -2 : .0001: 2;
 %Three frequencies to add together
freq1 = 2*pi*1; 
freq2 = 2*pi*2; 
freq3 = 2*pi*3;
 
 
sine1 = sin(freq1*time);
sine2 = sin(freq2*time); 
sine3 = sin(freq3*time);
sine4 =  sine1 + sine2 + sine3; 
 
cosine1 = cos(freq1*time);
cosine2 = cos(freq2*time);
cosine3 = cos(freq3*time);
cosine4 = cosine1 + cosine2 + cosine3;
 
figure(1)
subplot(321)
plot(time, sine1)
title('Sine wave of freq = 1 Hz')
axis([ -1 1 -3 3])
grid on
 
subplot(323)
plot(time, sine2)
title('Sine wave of freq = 2 Hz')
axis([ -1 1 -2 2])
grid on
 
subplot(325)
plot(time, sine3)
title('Sine wave of freq = 3 Hz')
axis([ -1 1 -2 2])
grid on
 
subplot(322)
plot(time, cosine1)
title('Cosine wave of freq = 1 Hz')
axis([ -1 1 -2 2])
grid on
subplot(324)
plot(time, cosine2)
title('Cosine wave of freq = 2 Hz')
axis([ -1 1 -2 2])
grid on
subplot(326)
plot(time, cosine3)
title('Cosine wave of freq = 3 Hz')
axis([ -1 1 -2 2])
grid on
 
figure(2)
subplot(311)
plot(time, sine4)
title('Sum of 3 Cosine waves')
axis([ -2 2 -3 4])
grid on
subplot(312)
 
plot(time, cosine4)
title('Sum of 3 Sine waves')
axis([ -2 2 -3 4])
grid on
 
subplot(313)
 
plot(time, cosine4+sine4)
title('Sum of both sine and cosines')
axis([ -2 2 -5 5])
grid on
 
%%
% Figure 1.10 – Sum of large number of sines and cosines together

N = 20; %number of sines or cosines to be summed
time = -2:.0001: 2;
sine_sum = 0;
cosine_sum = 0;

for n = 1:1: N
cosine_sum = cosine_sum + cos(2*pi*n*time)
end

subplot(311)
plot(time, (cosine_sum))
grid on
axis([ -2 2 -N/2 N*1.5])

title('(a) Sum of 20 harmonic cosine waves')
for n = 0:1: N;
sine_sum = sine_sum + sin(2*pi*n*time);
end

subplot(312)
plot(time, sine_sum)
grid on
axis([ -2 2 -N N])

title('(b) Sum of 20 harmonic sine waves')

subplot(313)
plot(time, (abs(cosine_sum))- (abs(sine_sum)))
grid on
axis([ -2 2 -20 N*1.5])
title('(c) Sum of 20 harmonic sine and cosine waves') % %% Odd and even square waves time = 0: .01: 5; f0 = 1*2*pi; %Manually sum the first four sinusiods to get odd and even square waves % Sine is an odd function so any waveform created from pure sines will also %be odd x1 = sin(f0 * time)+ (1/3)*sin(3*f0*time) + (1/5)*sin(5*f0*time)+ (1/7)*sin(7*f0*time); %Cosine is a even function so any waveform created from pure cosines will %be even x2 = cos(f0 * time)- (1/3)*cos(3*f0*time) + (1/5)*cos(5*f0*time)- (1/7)*cos(7*f0*time); subplot(211) plot(time, x1) title('Odd Square wave from Sines') axis([ 0 5 -1 1]) subplot(212) plot(time, x2) title('Even Square wave from Cosines') %% % Figure 1.15 Gibb's phenomenon % Gibbs time=linspace(-2,2,2000); square_wave=[zeros(1,500),2*ones(1,1000),zeros(1,500)]; k=2; N=[5,11,30,100]; for n=1:4; an=[]; for m=1:N(n) an=[an,2*k*sin(m*pi/2)/(m*pi)]; %calculate summed sine portion of square wave end; fN=k/2; for m=1:N(n) fN=fN+an(m)*cos(m*pi*time/2); %calculate summed cosine portion and add to sine end; title_string=int2str(N(n)); subplot(2,2,n),plot(time,square_wave,'--', 'LineWidth',.75); hold on; plot(time,fN,'-','LineWidth',1); hold off; axis([-2 2 -0.5 2.5]); grid; xlabel('Time'), ylabel('y_N(t)');title(['N= ',title_string]); end; % End of gibb.m %Chapter 2 %% %Figure 2.1 – Complex exponential %comexp produces the complex exponential diagram in Chapter 2. %Try changing function parameters to see the effect on the signal. clf time = 0:0.01:5; frequency=2*pi; %change this variable to see effect on complex exponential y=exp((j*(frequency))*time); figure(1); plot3(time,real(y),imag(y), 'linewidth', 0.3); grid xlabel('t'),ylabel('Re(y)'),zlabel('Im(y)'); %% % Figure 2.2 - Complex exponential with projections % Plot 3D complex exponential and also project the real and imaginary % parts onto their respective planes. % This will help you visualize what a complex exp. is made out of. clf time = 0:0.01:5; y1 = exp(1i*2*pi*time); y2 = exp(-1i*2*pi*time); subplot(121); plot3(time,real(y1), imag(y1),'b', 'linewidth', 1.5); %plot 3D complex exp grid xlabel('Time, t'),ylabel('Real plane'),zlabel('Imaginary plane'); th2 = linspace(-2,-2, length(real(y1))); %vector of -2 to hold projection on real plane hold on plot3(time, real(y1), th2, 'g', 'linewidth', 1.1); %plot real part against time, hold imaj part constant hold on th = linspace(2,2, length(real(y1))); %vector of 2 to hold projection on imaj plane plot3(time,th, imag(y1), 'r', 'linewidth', 1.1) %plot imaginary part of positive complex exp grid on % Negative complex exp subplot(122) plot3(time,real(y2), imag(y2), 'b', 'linewidth', 1.5); grid xlabel('Time, t'),ylabel('Real plane'),zlabel('Imaginary plane'); hold on plot3(time, real(y2), th2, 'g','linewidth', 1.1) hold on plot3(time,th, imag(y2), 'r', 'linewidth', 1.1) grid on % Chapter 3 %% % Figure 3.1 – Discrete sampling % Shows the digital vs. discrete signal % Discrete signals, better named as discrete time signals, can have sample points at any amplitude value. % This makes them discrete time but continuous amplitude. % Digital signals are limited in which amplitudes the sample can take. % Digital signals are discrete time and discrete amplitude. % Here we plot a continuous wave being sampled discretely and digitally clf time = 0: .01: 4; x = sin(2*pi*time)-.6*sin(2*pi*2.5*time)+.6; %continuous wave % setup the sampling frequency and period Fs = 8; ts = 1/Fs; n = 0: 32; xd = sin(2*pi*n*ts)-.6*sin(2*pi*2.5*n*ts)+.6; %discrete samples subplot(211) plot(time*Fs, x, '--') title('Exact Sampling') xlabel('Time') ylabel('Amplitude') axis([0 32 -1.5 2.5]) hold on stem(n, xd) subplot(212) plot(time*Fs, x, '--') title('Binary Digital Sampling') xlabel('Time') ylabel('Amplitude') axis([0 32 -1.5 2.5]) hold on % Set up the digital sampling by moving discrete time but continuous % amplitude samples to be discrete in both time and amplitude. for n=1:length(xd) if xd(n) > 0 xd(n) = 1 else xd(n) = -1 end end stem(linspace(0,32,33), (xd)) hold on %% % Figure 3.2 – Sampling a CT signal clf time = 0: .01: 2; x = .2*sin(2.5*pi*time)+ .9*cos(5*pi*time+1); Fs = 8; Ts = 1/Fs; k = 0: 19; xt = .2*sin(2.5*pi*k*Ts)+ .9*cos(5*pi*k*Ts+1); length(k) del = ones(length(k)); figure(1) plot(time*Fs, x, '-.') hold on xlim = get(gca,'xlim'); %Get x range hold on title('(a) - The analog signal') hFig = figure(1); % set(hFig, 'units', 'inches', 'pos', [7 5 4.5 1]) figure (2) stem(k, del) axis([ 0 20 0 1.25]) hFig = figure(2); % set(hFig, 'units', 'inches', 'pos', [7 5 4.5 1]) title('(a) - The Impulse train') figure(3) hFig = figure(3); % set(hFig, 'units', 'inches', 'pos', [7 5 4.5 1]) stem(k, xt) axis([ 0 20 -1.5 1.5]) hold on xlim = get(gca,'xlim'); %Get x range hold on plot([xlim(1) xlim(2)],[0 0],'k') plot(time*Fs, x, '-.') %% % square pulse response B = 16; tmin = -1/B; time = tmin*8: .01: -tmin*8; ht = sinc(B*time) figure(1) plot(time, (ht)); axis([ -8/B 8/B -.25 1.25]) xlabel('Time, seconds') grid on title('(a) - Pulse applied to every sample, width =.075 sec.') %Frequency domain figure(2) f = linspace(2/tmin, -2/tmin, length(time)); y = (1/B)*rectpuls(f, B) plot(f, y) title('(b) - Frequency response of the sinc pulse, B = 2/T = 8 Hz') xlabel('Frequency, Hz') grid on %% % Fig. 3.5 - Zero and first-order hold % Reconstructs a sampled signal with two methods. % Zero order hold will maintain the same amplitude as previous sample until % next sample point. % First order will linearly interpolate from one sample to the next. clf f=4; %frequency of the impulse in Hz fs=8; % sample frequency is 10 times higher time=-4:1/fs:3; % discrete time vector y2 = ones(size(time)); plot(time,y2); subplot(211) time2 = -4: .01:3; % "continuous" time vector xc_discrete = cos(2*pi*time+pi/3)*.2 - sin(2*pi*2.*time+pi/3)*.8 - .01*time; xc2 = cos(2*pi*time2+pi/3)*.2 - sin(2*pi*2.*time2+pi/3)*.8 - .01*time2; stairs(time, xc_discrete, 'r') % creates the zero order hold hold on plot(time2, xc2, '-.') title('Zero Order Reconstruction') xlabel('Time') ylabel('Amplitude') grid on axis([ -1 1.6 -1 1]) set(gca, 'Xtick', -1: .125:1.6) subplot(212) plot(time, xc_discrete, 'r') hold on plot(time2, xc2, '-.') %Matlab's plot function plots discrete samples by using a first order hold! title('First Order Reconstruction') xlabel('Time') ylabel('Amplitude') grid on axis([ -1 1.6 -1 1]) set(gca, 'Xtick', -1: .125:1.6) %% % Fig. 3.6 - Sinc reconstruction % clf f=4; %frequency of the impulse in Hz fs=8; % sample frequency is 10 times higher time=-4:1/fs:3; % time vector y2 = ones(size(time)); plot(time,y2); subplot(311) n = 0: 8 time2 = -4: .01:3; xc = cos(2*pi*time+pi/3)*.2 - sin(2*pi*2.*time+pi/3)*.8 - .01*time; xc2 = cos(2*pi*time2+pi/3)*.2 - sin(2*pi*2.*time2+pi/3)*.8 - .01*time2; xs = y2.*xc; plot(time, xs, 'o') hold on plot(time2, xc2, '-.') grid on axis([ -1 1.6 -1 1]) set(gca, 'Xtick', -1: .125: 1.6) subplot(312) n2 = 0:56; T = 1/fs; fb = 1; h = 0; for i = 1: length(xs) t3 = -8: .01: 32; h = xs(i)*sinc((t3-n2(i))); plot(t3-8, h) hold on end axis([ -8 13 -1 1]) set(gca, 'Xtick', -8: 1: 13) grid on hold on plot(time*8, xs, 'o') subplot(313) h = 0; for i = 1: length(xs) t3 = -8: .01: 32; h = h + xs(i)*sinc((t3-n2(i))); hold on end plot(t3-8, h) axis([ -8 13 -1 1]) grid on %% % Figure 3.9 – Sampled signals % This plots a continous wave with continuous time. % Underneath, the sampled waveform is plotted with X-axis in terms of % sample number and in terms of radians. For a discretely sampled signal, % the sampling rate and the sampled singal's fundamental frequency will % determine how much each sample is in terms of radians. cls clf time = -.5: .01: .5; %continuous time, to be convertered to digital freq and samples y1 = cos(8*pi*time); %continuous cosine signal to be sampled clf; subplot(3,1,1) plot(time, y1, 'linewidth', 1.00) grid on; xlabel('Time, t seconds'); ylabel('(a) - x(t)'); axis; subplot(3,1,2) n = -5: 5; %sample number vector y2 = cos(4*pi*n*.2); % sampled version of continuous cosine stem(n, y2, 'filled') set(gca,'xTick',-5: 1: 5) %axis here is sample number grid; xlabel('Sample, k'); ylabel('(b) - x[k]'); subplot(3,1,3) n = -5: 1: 5; n2 = -2*pi: 2*pi/5: 2*pi; y2 = cos(2*n2); stem(n2, y2, 'filled') axis([-2*pi 2*pi -1 1]); grid; xlabel('Radians'); ylabel('(c) - x[k]'); axis([-2*pi 2*pi -1. 1.]) % Define x-ticks and their labels to be in terms of digital freq set(gca,'xTick',-2*pi: pi/5: 2*pi) set(gca,'xTickLabel',{'-2pi', '', '-8pi/5', '', '-6pi/5', '', '-4pi/5', '', '-2pi/5', '', '0', '', '2pi/5', '', '4pi/5', '', '6pi/5', '', '8pi/5' , '', '2pi'}) %% % Figure 3.11 three signals with same samples. % We plot 3 continuous waves that when sampled at a given rate, end up % having the same samples. The signals are aliases at this sampling rate. clf freq1 = 2 freq2 = 6; freq3 = 10; time = 0: .001: 1; x1 = cos(2*pi*freq1*time); x2 = cos(2*pi*freq2*time); x3 = cos(2*pi*freq3*time); Fs = 4; n = 0: 8; x4 = cos(2*pi*n/Fs); plot(time, x1,'r') hold on plot(time, x2,'b') plot(time, x3,'g') stem(n/(2*Fs), x4, 'filled', 'linewidth', .2) axis([ 0 1 -1 1]) xlim = get(gca,'xlim'); %Get x range hold on plot([xlim(1) xlim(2)],[0 0],'k') title('Three different waves, f = 2, 6, 10 Hz all pass through same samples') xlabel('Time, t') %% % Figure 3.6 – Sinc reconstruction % A waveform is sampled to bring it to discrete time domain. It can be brought back by % placing a scaled sinc at each sample point and summing all sincs together. clf % set(0,'DefaultLineMarkerSize', 3) % set(0,'DefaultLineLineWidth', 1.0) % f = 1: .01: 1; time = 0: .01: 1; %waveform to be sampled and reconstructed x = .3*sin(2*pi*time)-.6*cos(6*pi*time); Fs = 10; k = 0: 30; %discrete samples of waveform xn = .3*sin(2*pi*k/Fs)-.6*cos(6*pi*k/Fs); figure(1) plot(time, x) xlabel('Time') ylabel('Amplitude') title('A sampled signal') hold on stem(k/Fs, xn, 'filled') axis([ 0 1 -1 1]) figure(2) clf % Create sinc for each sample point %Shift each sinc by n/10 and scale it by sample amplitude ht1 = xn(1)*sinc(Fs*((time-(0/10)))); ht2 = xn(2)*sinc(Fs*((time-(1/10)))); ht3 = xn(3)*sinc(Fs*((time-(2/10)))); ht4 = xn(4)*sinc(Fs*((time-(3/10)))); ht5 = xn(5)*sinc(Fs*((time-(4/10)))); ht6 = xn(6)*sinc(Fs*((time-(5/10)))); ht7 = xn(7)*sinc(Fs*((time-(6/10)))); ht8 = xn(8)*sinc(Fs*((time-(7/10)))); ht9 = xn(9)*sinc(Fs*((time-(8/10)))); ht10 = xn(10)*sinc(Fs*((time-(9/10)))); ht11 = xn(11)*sinc(Fs*((time-(10/10)))); %plot all sincs plot(time, ht1, time, ht2, time, ht3, time, ht4, time, ht5, time, ht6, time, ht7, time, ht8, time, ht9, time, ht10, time, ht11) hold on stem(time(1: 10: 100), xn(1:1:10), 'filled') hold on xlabel('Time') ylabel('Amplitude') title('Each sample in one period replaced by scaled sinc function') xlim = get(gca,'xlim'); %Get x range plot([xlim(1) xlim(2)],[0 0],'k', 'linewidth', .75) %Plot sinc reconstruction and original signal figure(3) hta = ht1+ht2+ht3+ht4+ht5+ht6+ht7+ht8+ht9+ht10+ht11; plot(time, hta, '--') hold on xlim = get(gca,'xlim'); %Get x range plot([xlim(1) xlim(2)],[0 0],'k', 'linewidth', .75) hold on plot(time, x,'r') xlabel('Time') title('Summation of all scaled sinc functions vs. original signal') axis( [ 0 1 -1.5 1.5]) time = -1: .01: 1; %Show sinc basis functions on their own figure(4) ht0 = xn(1)*(Fs/pi)*sinc(Fs*time/pi); mht = max(-ht0); ht0=(-1)* ht0/mht; plot(time, ht0, time+1/10, ht0, time+2/10, ht0, time+3/10,ht0, time+4/10, ht0, time+5/10, ht0) xlabel('Time') title('Sinc functions ready to replace samples') hold on xlim = get(gca,'xlim'); %Get x range plot([xlim(1) xlim(2)],[0 0],'k', 'linewidth', .75) %% % Example 3.7 clf k = -20: 19 Fs = 10; subplot(3,1,1) x = 1+sin(2*pi*k/Fs); stem(k, x, 'filled') xlabel('Sample, k') grid on axis([ -20 20 0 3]) ffx = fftshift((1/10)*abs(fft(x, 10))); n = -5: 4; subplot(3,1,2) stem(n, ffx, 'filled') axis([ -5 5 0 1.5]) xlabel('Bin') set(gca,'xTick',-5: 5) grid on newx = [ ffx ffx ffx ffx ffx] time2 = linspace(-25, 24, 50) subplot(3,1,3) stem(time2, newx, 'filled') xlabel('Bin') grid on axis([ -25 25 0 1.5]) %% % Figure 3.22 % Here we see single and replicating spectrums clf n = -10: 9; % sample number Fs = 10; % sampling frequency figure(1) %x = [0 1 2 3]; x = 1+sin(2*pi*n/10); k2 = -length(n)*2.5: length(n)*2.5-1; x1 = [ x x x x x] figure(1) stem(k2, x1, 'filled') xlabel('Sample, n') ylabel('x[n]') title('(a) - the discrete signal') axis([ -10 10 -1 3]) ffx = fftshift((1/10)*abs(fft(x, 10))) n = -5: 4; figure(2) stem(n, ffx, 'filled') axis([ -5 4 0 1.5]) xlabel('Frequency') ylabel('C[k]') title('(b) - Fundamental Spectrum') grid on figure(3) newx = [ ffx ffx ffx ffx ffx ] time2 = linspace(-25, 24, 50) stem(time2, newx, 'filled') xlabel('Frequency') ylabel('C[k]') title('(c) - Repeating coefficients') axis([ -25 25 0 1.5]) %% % Figure 3.23 % When calculating the discrete fourier series, there are an infinite number of frequencies % that can go through the sample points. These frequencies are called % aliases. They show up as repeating bands reproduced at the sampling % frequency. Here we plot the center band, called the fundamental spectrum % and a wider view showing more alias bands. clf n = -20: 19 % sample index Fs = 10; % sampling freq - this will set the spacing of alias spectrums figure(1) x = .5 +.25*cos(2*pi*n/5) - .6*sin(pi*n/2); n2 = -length(n)*1.5: length(n)*1.5-1; set(0,'DefaultLineMarkerSize', 3) x1 = [ x x x]; s(1)=subplot(3,1,1); stem(n2, x1, 'filled') xlabel('Sample, n') axis([ -20 20 0 2]) ffx = fftshift((1/4)*abs(fft(x, 20))); n = -10: 9; grid on s(2)=subplot(3,1,2); stem(n, ffx, 'filled') axis([ -10 10 0 3]) xlabel('Bin') grid on s(3)=subplot(3,1,3); newx = [ ffx ffx ffx ffx ffx]; time2 = linspace(-50, 49, 100); stem(time2, newx, 'filled') xlabel('Bin') axis([ -50 50 0 3]) grid on set(gca,'xTick',-50:10: 50) title(s(1),'Time domain samples') title(s(2),'Frequency domain, fundamental alias') title(s(3),'Frequency domain, multiple alias spectrums') %% %Figure 3.21 – Real and imaginary components of the discrete signal harmonics. % This figure plots the real parts of a discrete sampled exponential %Discrete harmonics repeat every 2*pi while continuous harmonics do not repeat. %Note that while the discrete harmonic frequency increases the number of %samples stays the same. n=-6:6; %time axis set(0,'DefaultLineMarkerSize', 3) set(0,'DefaultLineLineWidth', .75) clf N=6; %number of harmonics to divide the period by w0=2*pi/N; %digital frequency axis([-12.5 12.5 -1.1 1.1]); % 0th Harmonic k=0; Phi0n=exp(j*w0*k*n); subplot(3,3,1) stem(n,linspace(1,1,13), 'filled'); hold on plot(time, cos(w0*k*time)) axis([-6 6 -1.1 1.1]); title('2pi(0/6)') hold on time= -1: 1/(2*N): 1; plot(time*N, cos(w0*k*time)) axis([-6 6 -1.1 1.1]); k=1; Phi1n=exp(j*w0*k*n); subplot(3,3,2); stem(n,real(Phi1n), 'filled'); hold on title('2pi(1/6)') time= -1*N: 1/(2*N): 1*N; plot(time, cos(w0*k*time)) axis([-6 6 -1.1 1.1]); k=2; Phi2n=exp(j*w0*k*n); subplot(3,3,3); stem(n,real(Phi2n), 'filled'); hold on title('2pi(2/6)') time= -1*N: 1/(2*N): 1*N; plot(time, cos(w0*k*time)) axis([-6 6 -1.1 1.1]); k=3; Phi3n=exp(j*w0*k*n); subplot(3,3,4); stem(n,real(Phi3n), 'filled'); title('2pi(3/6)') hold on time= -1*N: 1/(2*N): 1*N; plot(time, cos(w0*k*time)) axis([-6 6 -1.1 1.1]); k=4; Phi4n=exp(j*w0*k*n); subplot(3,3,5);stem(n,real(Phi4n), 'filled'); title('2pi(4/6)') hold on time= -1*N: 1/(2*N): 1*N; plot(time, cos(w0*k*time)) axis([-6 6 -1.1 1.1]); k=5; Phi5n=exp(j*w0*k*n); subplot(3,3,6); stem(n,real(Phi5n), 'filled'); title('2pi(5/6)') hold on time= -1*N: 1/(2*N): 1*N; plot(time, cos(w0*k*time)) axis([-6 6 -1.1 1.1]); %Real part of 6th Harmonic will be repeat of 0th in discrete domain k=6; Phi6n=exp(j*w0*k*n); subplot(3,3,7); stem(n,real(Phi6n), 'filled'); title('2pi(6/6)') hold on time= -1*N: 1/(2*N): 1*N; plot(time, cos(w0*k*time)) axis([-6 6 -1.1 1.1]); %Real part of 7th Harmonic will be repeat of 1st in discrete domain k=7; subplot(3,3,8); Phi7n=exp(j*w0*k*n); stem(n,real(Phi7n), 'filled'); axis([-6 6 -1.1 1.1]); hold on plot(time, cos(w0*k*time)) title('2pi(7/6)') %Real part of 8th Harmonic will be repeat of 2nd in discrete domain k=8; subplot(3,3,9); Phi8n=exp(j*w0*k*n); subplot(3,3,9); stem(n,real(Phi8n), 'filled'); hold on plot(time, cos(w0*k*time)) title('2pi(8/6)') axis([-6 6 -1.1 1.1]); %% %Figure 3.22 close all k = -20: 19 Fs = 10; subplot(3,1,1) x = 1+sin(2*pi*k/Fs); stem(k, x, 'filled') xlabel('Sample, k') grid on<br>axis([ -20 20 0 3]) ffx = fftshift((1/10)*abs(fft(x, 10))); n = -5: 4; subplot(3,1,2) stem(n, ffx, 'filled') axis([ -5 5 0 1.5]) xlabel('Bin') set(gca,'xTick',-5: 5) grid on newx = [ ffx ffx ffx ffx ffx]<br>t2 = linspace(-25, 24, 50) subplot(3,1,3) stem(t2, newx, 'filled') xlabel('Bin') grid on %% % Figure 3.24 close all k = -10: 9 Fs = 10; subplot(3,1,1) x = [0 0 2 1]; x1 = [ x x x x x] stem(k, x1, 'filled') xlabel('Sample, k') axis([ -10 10 0 3]) ffx = fftshift((1/4)*abs(fft(x, 4))) n = -2: 1; grid on subplot(3,1,2) stem(n, ffx, 'filled') axis([ -3 3 0 1.0]) xlabel('Bin') grid on subplot(3,1,3) newx = [ ffx ffx ffx ffx ffx] t2 = linspace(-10, 9, 20); stem(t2, newx, 'filled') xlabel('Bin') grid on axis([ -10 10 0 1.0]) %% % Figure 3.25 clf k = -10: 9 Fs = 10; figure(1) x = [0 1 2 3]; x1 = [ x x x x x] subplot(3,1,1) stem(k, x1, 'filled') xlabel('Sample, n') grid on axis([ -10 10 0 3]) ffx = fftshift((1/4)*abs(fft(x, 4))) n = -2: 1; set(gca,'xTick',-10: 4: 10) subplot(3,1,2) stem(n, ffx, 'filled') axis([ -3 3 0 2]) xlabel('Bin') grid on set(gca,'xTick',-3: 3) subplot(3,1,3) newx = [ ffx ffx ffx ffx ffx] time2 = linspace(-10, 9, 20) stem(time2, newx, 'filled') xlabel('Bin') grid on axis([ -10 10 0 2]) set(gca,'xTick', -10 : 4: 10) %% % Figure 3.26 - Square pulse example fig=gcf; % set(findall(fig,'-property','FontSize'),'FontSize',8) % % Change default axes fonts. % set(0,'DefaultAxesFontName', 'charter') % set(0,'DefaultAxesFontSize', 8) % set(0, 'DefaultAxesColorOrder', [ 0 0 0]) % set(0, 'DefaultAxesLineStyleOrder', '-|:|--|-.') % set(0,'DefaultAxesLineWidth', .75) % set(0,'DefaultLineLineWidth', 1.0) % set(0,'DefaultLineMarkerSize', 2) % set(0,'DefaultTextFontname', 'charter') % set(0,'DefaultTextFontSize', 8) clf L = 2; N = 10; k = [-N/2: N/2]; n = [0:1:N-1]; k=[0:1:N-1]; xn = [ones(1,L), zeros(1,N-L)]; xn2 = [ xn xn xn xn]; subplot(3,1,1) lk = linspace(0, 4*N-1, 4*N) stem(lk, xn2,'filled'); %str=sprintf('Square pulse with N = %d , L = %d', N, L); %title(str) subplot(3,1,2) WN = exp(-j*2*pi/N); nk = n'*k; k3 = linspace(-5,4, 10) WNnk = WN.^nk; Xk = xn * WNnk; stem(k3, fftshift(abs(Xk)),'filled') hold on plot(k3, fftshift(abs(Xk)), '--') axis([ -5 , 4, -.5, N]) subplot(3,1,3) lk4 = linspace(-20, 19, 40) Xn2 = [ Xk Xk Xk Xk]; stem(lk4, Xn2, 'filled'); hold on plot(lk4, Xn2, '--') axis([-20, 19, -1, max(Xn2)]) % Chapter 4 %% % sampling a signal Close all time = 0: .01: 45 y = sin((1/14)*2*pi*time)+ .3*sin((1/14)*2*pi*time*3) ; figure(1) plot(time, y, 'linewidth', 1.5) xlabel('Time, t') ylabel ('Amplitude') title('Continuous signal') grid on hold on plot([xlim(1) xlim(2)],[0 0],'k', 'linewidth', .75) ts = 1/14; n = 0: 1: 44; yd = sin(2*pi*n*ts) + .3*sin(2*pi*ts*n*3); stem(n, yd, 'filled') xlabel('Sample, k') ylabel ('Amplitude') title('Sampled signal') % Chapter 5 %% % Figure 5.6 - DTFT of a group of pulses clf Kp = 7; b = [0 0 0 0 ones(Kp, 1)' 0 0 0 0 ]; lenb = length(b); subplot(121) n2 = -4-(Kp-1)/2: 4+(Kp-1)/2; stem(n2, b, 'filled') grid on axis([-4-(Kp-1)/2 4+(Kp-1)/2 0 1.5]) subplot(122) time = -10:.01: 10; xom = abs(diric(2*time, Kp)); grid on plot(time, xom) axis([-3*pi 3*pi -2 2]) grid on set(gca,'XTick',[-3*pi:pi/2:3*pi]) pilabels %% %% Figure 5.7 – Sinc and Drichlet functions clf N = 10; k = -N/2: N/2; om = -N/2: .1: N/2; x = sinc(k); x1= sinc(om); subplot(221) stem(k, x) axis([-N/2, N/2, -.5 1.5]) set(gca,'XTick',[-N/2:N/2]) grid on subplot(223) stem(k, x) hold on plot(om, (x1), '-.') grid on axis([-N/2, N/2, -.5 1.5]) set(gca,'XTick',[-N/2:N/2]) subplot(224) stem(k, x) hold on plot(om, abs(x1), '-.') grid on axis([-N/2, N/2, -.5 1.5]) set(gca,'XTick',[-N/2:N/2]) subplot(222) N = 50; k = -N/2: N/2; om = -N/2: .1: N/2; x = sinc(k); stem(k, x) axis([-N/2, N/2, -.5 1.5]) grid on %% %Figure 5.6 close all Kp = 7; b = [0 0 0 0 ones(Kp, 1)' 0 0 0 0 ]; lenb = length(b); subplot(121) n2 = -4-(Kp-1)/2: 4+(Kp-1)/2; stem(n2, b, 'filled') grid on axis([-4-(Kp-1)/2 4+(Kp-1)/2 0 1.5]) subplot(122) t = -10:.01: 10; xom = abs(diric(2*t, Kp)); grid on plot(t, xom) axis([-3*pi 3*pi -2 2]) grid on set(gca,'XTick',[-3*pi:pi/2:3*pi]) pilabels %% % Figure 5.7 - The Drichlet function clf cls x = linspace(-3*pi,3*pi,300); subplot(221); plot(x,(diric(x,2))); title('Diric, N = 2'); axis tight; pilabels grid on ylabel('Amplitude') subplot(222); plot(x,(diric(x,3))); title('Diric, N = 3'); axis tight; pilabels grid on subplot(223); plot(x,abs(diric(x,5))); title('Diric, N = 5'); axis tight; pilabels ylabel('|Amplitude|') xlabel('Frequency') grid on subplot(224); plot(x,abs(diric(x,8))); title('Diric, N = 8'); axis tight; pilabels grid on xlabel('Frequency') %% % Figure 5.8 cls Clf n2 = -30: 30; time = -40: .01: 40; T0 = 1; tau = .25; xfscn = [ 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 1]; xfsc = sinc(tau*time/T0); xfscn = fft(xfscn, 20) n3 = -10: 1: 9; plot(n3, fftshift(abs(xfscn))) figure plot(n3, imag(xfscn)) %% clf K0 = 8; N = 4; a = [1 ] b = [ a 0 0 0 0 0 0 ]; c = [ b b b]; subplot(2,1,1) time = -11: 9; stem(time, c) axis([-11 9 0 1.5]) subplot(2,1,2) n2 = -14: 13; xom = (N/K0)*diric(2*n2*pi/K0, N); grid on stem(n2, xom, 'filled') hold on n2 = -11: .01: 9; xom2 = (N/K0)*diric(2*n2*pi/K0, N); plot(n2, xom2, '-.r' ) axis([-11 9 -1 1]) %% % CTFT %Compare all % Here we plot a time domain signal, then its CTFT, its phase, and the % real and imaginary parts. clf % Create a "continuous" rectangle pulse signal % Plot it x = [zeros(1, 124) ones(1,8) zeros(1,124)]; time = -128: 127; figure(1) plot(time, x) grid on axis([ -128 128 0 1.5]) % On this figure we plot the CTFT and its parts (phase, real, imag) figure(2) subplot(2,2,1) % We use the FFT to plot the CTFT. Technically it is the DTFT but with % enough sample points, the two look the same when plotted only between % -pi/2 to pi/2 ctft = fft(x, 256) % FFTshift rearranges the array of FFT values to be centered at 0 going to % +/- (pi/2). Makes visualizing the FFT much simpler. x2 = fftshift(abs(ctft))' % Create the index by which we plot the CTFT against. freq = linspace(-pi/2,pi/2,length(time)); plot(freq, x2) axis([ -pi/2 pi/2 0 10]) % set(gca,'XTick',[-128:64:128]) title('Abs of ctft') grid on subplot(2,2,2) plot(freq, (imag(ctft))); % set(gca,'XTick',[-128:64:128]) %set(gca,'XTickLabel',['-128'; '-64 '; '0'; ' 64'; '128']) phase = atan(imag(ctft)./(real(ctft))) plot(freq, phase) axis([ -pi/2 pi/2 -pi pi]) set(gca,'YTick',[-pi:pi/4:pi]) title('Phase') grid on set(gca,'XTick',[-128:64:128]) %set(gca,'XTickLabel',['-128'; '-64 '; '0'; ' 64'; '128']) subplot(2,2,3) plot(time, real(fftshift(ctft))) title('Real part') axis([ -128 128 -10 10]) grid on set(gca,'XTick',[-128:64:128]) %set(gca,'XTickLabel',['-128'; '-64 '; '0'; ' 64'; '128']) subplot(2,2,4) plot(time, imag(ctft)) title('Imaginary part') axis([ -128 128 -1.5 1.5]) set(gca,'XTick',[-128:64:128]) %set(gca,'XTickLabel',['-128'; '-64 '; '0'; ' 64'; '128']) grid on %% time = -24: 24; x = time.^2/10; plot(time, x) set(gca,'XTick',[-24:8:24]) set(gca,'XTickLabel',['-24'; '-16 '; '-8'; ' 0'; '8'; '16 '; '24']) %% %% Figure 5.5 – DTFT of a discrete signal clf x = [0 2 -3 0 5 zeros(12, 1)' ]; n = 0: 16; figure(1) stem(n, x) grid on axis([ 0 16 -3 6]) f = -1: .01: 1; ce = 2*exp(-1i*2*pi*f*1)-3*exp(-1i*2*pi*f*3)+5*exp(-1i*2*pi*f*5) xf = abs(ce); figure(2) plot(f, xf) grid on figure(3) phase = atan((imag(ce)./real(ce))) plot(f, phase) grid on figure(3) k = 3; cen = 0; for n = 0:4 cen = abs(cen) + 2*exp(-1i*2*pi*k*n)-3*exp(-1i*2*pi*k*n)+5*exp(-1i*2*pi*k*n) end %% % Diric function clf N = 10; NT = 48; k = -N/2: N/2; x = abs(sinc(k)) figure(1) grid on %stem(k, x) hold on om = -4*pi: 2*pi/NT: 4*pi %plot(om, abs(sinc(om))) axis([ -2*pi 2*pi 0 1]) grid on subplot(322) xd = abs(diric(om, 3)) plot(om, xd) axis([ -3*pi 3*pi 0 1]) subplot(321) plot(om, xd) axis([ -pi pi 0 1]) text(-2.5, 0.8, 'N = 3') subplot(324) xd2 = abs(diric(om, 5)) plot(om, xd2) axis([ -3*pi 3*pi 0 1]) subplot(323) plot(om, xd2) axis([ -pi pi 0 1]) text(-2.5, 0.8, 'N = 5') subplot(326) xd2 = abs(diric(om, 9)) plot(om, xd3) axis([ -3*pi 3*pi 0 1]) subplot(325) plot(om, xd3) axis([ -pi pi 0 1]) text(-2.5, 0.8, 'N = 9') %% % Figure 5.11 Rs = 1; %Sets the frequency of source. Ts = 1/Rs; %Calculates the period of the source. al = .5; f = 0:0.01:T; %Defines vector t with 10T+1 elements from 0-T. s = zeros(size(f)); for i = 1:length(f) if f(i)<(1-al)/(2*Ts) s(i) = Ts; elseif (1-al)/(2*Ts)<=f(i)& f(i)<(1+al)/(2*Ts) s(i) = (Ts*.5)*(1+cos((pi*Ts/al)*(abs(f(i))-((1-al)/(2*Ts))))); else s(i) = 0; end end figure(1) %Names the figure plot(f,s); %Draws a plot where t is the x-axis and s(t) is the y-axis xlabel('time (seconds)'); %Inserts a x-axis label on the plot ylabel('s(t)'); %Inserts a label for the y-axis grid on hold on %% % Triangualar pulse x = [ zeros(7,1)' .25 .5 .75 1 .75 .5 .25 zeros(6,1)']; lenx = length(x) n = -lenx/2: lenx/2-1; subplot(311) stem(n, x) grid on N = 24; om = 2*pi/N; xf = fftshift(abs(fft(x, 20))/20); n2 = -10: 9; subplot(312) plot(n2, xf) grid on xfr = [ xf xf xf xf xf] subplot(313) n3 = -50: 49 n4 = -5*pi: pi/10: 5*pi-pi/10; plot(n4, xfr) grid on set(gca,'XTick',[-10*pi:2*pi/2:10*pi]) pilabels %% % Figure 5.10 – Time expansion property a = [ 0 2 -3 0 5 0]; x = [ a zeros(12, 1)']; lx = length(x) n = -lx/2 : lx/2-1; subplot(211) stem(n,x) grid on subplot(212) om = -3*pi: .01: 3*pi; xf = 2*exp(j*om)-3*exp(j*2*om)+5*exp(j*4*om); plot(om, abs(xf)) grid on pilabels %% % Example time expansion a = [ 1 0 0 1 0 0 1]; x = [ zeros(5, 1)' a zeros(4, 1)']; lx = length(x) n = -lx/2 : lx/2-1; subplot(221) stem(n,x) grid on subplot(222) om = -6*pi: pi/6: 6*pi; xf = diric(om/3, 3); plot(om, abs(xf)) set(gca,'XTick',[-6*pi:1*pi:6*pi]) grid on pilabels a = [ 1 1 1]; x = [ zeros(6, 1)' a zeros(7, 1)']; lx = length(x) n = -lx/2 : lx/2-1; subplot(223) stem(n,x) grid on subplot(224) om = -6*pi: pi/6: 6*pi; xf = diric(om/1, 3); plot(om, abs(xf)) grid on pilabels set(gca,'XTick',[-6*pi:1*pi:6*pi]) %% Sine om0 = 1; n2 = -10: 10; x = sin(om0*n2) lx = length(x) n = -lx/2 : lx/2-1; subplot(211) stem(n,x) axis([ -10 10 -2 2]) grid on subplot(212) om = -3*pi: .01: 3*pi; NF = 64; xf = fftshift(fft(x, 64)*1/64) n3 = -NF/2: NF/2-1 stem(n3, real(xf)) grid on axis([ -NF/2: NF/2-1 -1 2]) %% % Figure 5.13 - Gaussian pulse DTFT close all n = -20: 20; sig = 4; x = (1/(sig*2*pi)^.5)*exp(-n.^2/(2*sig)); subplot(211) stem(n, x) grid on subplot(212) om = -4*pi:pi/8: 4*pi xf = exp(-om.^2/(2*sig)); plot(om, xf) grid on pilabels % this can be found on Matlab file Exchange %% ND = 11; om = -1*pi: .01: 1*pi; xd3 = 10*log10(abs(diric(om, ND)).^2) plot(om, xd3) set(gca, 'Xtick', [-pi: 2*pi/ND: pi]) grid on %pilabels for k = 1: 5 x1 = diric((k * pi/ND), ND); x1 = 20*log10(abs(x1)) end %% %Examples of DTFT om = -3*pi : .01 : 3*pi; n = -6: 6 x1 = [ 0 0 0 0 0 1 1 1 0 0 0 0 0]; x2 = [ 0 0 0 0 1 1 1 1 1 0 0 0 0]; subplot(221) stem(n, x1); ylabel('Amplitude') title('(a) x[n] = 111') grid on subplot(222) stem(n, x2); title('(b) x[n] = 11111') grid on x1 = (exp(1i*om*1)) + 1 + (exp(-1i*om*1)); x2 = 1 + exp(1i*om*2)+ exp(1i*om*1)+ exp(-1i*om*1)+ exp(-1i*om*2); subplot(223) plot(om, (x1)) title('(c) DTFT - 3 impulses') grid on pilabels xlabel('Frequency') ylabel('Amplitude') subplot(224) plot(om, (x2)) title('(d) DTFT - 5 impulses') grid on pilabels xlabel('Frequency') %% %% b = [ 42 39 36 33 30 27 24 21 18 15]; CI = [ 1 1 4 13 19 31 50 62 18 2]; CI = CI./201; bar(b, CI) xlabel('CI') ylabel(' Frequency of occurence') %% % Figure 5.16 – DTFT of a periodic signal clf NT = 64; k = -N/2: N/2; om = -4*pi: .01: 4*pi subplot(321) ND = 5; om = -4*pi: .01: 4*pi xd3 = (diric(om, ND)) plot(om, xd3) axis([ -pi pi -.5 1]) hold on f = -4*pi : 2*pi/ND: 4*pi xdi = diric(f, ND) stem(f, xdi) axis([ -pi pi -.5 1]) subplot(322) ND = 5; om = -4*pi: .01: 4*pi xd3 = (diric(om, ND)) plot(om, xd3) axis([ -12 12 -.5 1]) hold on f = -4*pi : 2*pi/ND: 4*pi xdi = diric(f, ND) stem(f, xdi) subplot(323) ND = 7; om = -1*pi: .01: 1*pi xd3 = (diric(om, ND)) plot(om, xd3) axis([ -pi pi -.5 1]) hold on f = -4*pi : 2*pi/ND: 4*pi xdi = diric(f, ND) stem(f, xdi) subplot(324) ND = 7; om = -4*pi: .01: 4*pi xd3 = (diric(om, ND)) plot(om, xd3) axis([ -12 12 -.5 1]) hold on f = -4*pi : 2*pi/ND: 4*pi xdi = diric(f, ND) stem(f, xdi) subplot(325) ND = 9; om = -1*pi: .01: 1*pi xd3 = (diric(om, ND)) plot(om, xd3) axis([ -pi pi -.5 1]) hold on f = -4*pi : 2*pi/ND: 4*pi xdi = diric(f, ND) stem(f, xdi) subplot(326) ND = 9; om = -4*pi: .01: 4*pi xd3 = (diric(om, ND)) plot(om, xd3) axis([ -12 12 -.5 1]) hold on f = -4*pi : 2*pi/ND: 4*pi xdi = diric(f, ND) stem(f, xdi) %% example 2 a = [ 1 2 4 0]; x = [ a zeros(12, 1)']; lx = length(x) n = -lx/2 : lx/2-1; subplot(211) stem(n,x) grid on subplot(212) om = -3*pi: .01: 3*pi; xf = 1+ 2*cos(om*2)+4*cos(4*om)-j*(sin(2*om)+2*sin(4*om)); plot(om, abs(xf)) grid on % Chapter 6 %% % Example 6.1 Simple DTFT clf beep off x0 = [ 0 2 -3 0 5 zeros(12, 1)' ]; x = [ 0 2 -3 0 5 ]; lx = length(x); n = 0: lx-1; subplot(311) stem(n, x) grid on axis([ 0 lx -4 6]) om = -pi: .01: pi; ce = 2*exp(-1i*om*1)-3*exp(-1i*om*2)+5*exp(-1i*om*4) xf = abs(ce); grid on subplot(312) N = lx; xdft = fftshift(abs(fft(x, N))) Nn = (-lx+1)/2: (lx-1)/2 stem(Nn, xdft) grid on subplot(313) stem(Nn*(2*pi/lx), xdft) hold on grid on plot(om, xf, '-.') %% Chapter 7 %% Periodgram for small N and large with windowing freq1 = 64; n = -1024: 1023; clf n0 = 0.05*randn(length(n),1); fs = 8; %data = xlsread('IQdata',1, 'a1:A2048'); % QPSK signal out of non-linear amplifier n = 1: 2048; n0 = + .01*randn(2048, 1)'; data = sin(2*pi*n/fs)+.6*cos(2*pi*2.4*n/fs); subplot(221) NFFT = 64; x = data(1:NFFT); w = 1; x = w.* x'; xf = fft(x,NFFT); Pxx=10*log10(xf.*conj(xf)); Pxx = Pxx - max(Pxx); f = (-NFFT/2)*fs/NFFT: fs/NFFT: (NFFT/2-1)*fs/NFFT; %Frequency Vector plot(f,fftshift(Pxx), 'linewidth', .25); axis([ -4 4 -80 0]) title('(a) Periodogram, No window ') ylabel('PSD, dB') set(gca, 'Xtick', -4: 2: 4) grid on subplot(222) NFFT = 64; x = data(1:NFFT); w = hanning(length(x(1:NFFT))) x = w.* x'; xf = fft(x,NFFT); Pxx=10*log10(xf.*conj(xf)); Pxx = Pxx - max(Pxx); f = (-NFFT/2)*fs/NFFT: fs/NFFT: (NFFT/2-1)*fs/NFFT %Frequency Vector plot(f,fftshift(Pxx), 'linewidth', .5); axis([ -4 4 -80 0]) title('(b) Periodogram, Hanning') ylabel('PSD, dB') set(gca, 'Xtick', -4: 2: 4) grid on subplot(223) NFFT = 64; x = data(1:NFFT); w = blackman(length(x(1:NFFT))) x = w.* x'; xf = fft(x,NFFT); Pxx=10*log10(xf.*conj(xf)); Pxx = Pxx - max(Pxx); f = (-NFFT/2)*fs/NFFT: fs/NFFT: (NFFT/2-1)*fs/NFFT; %Frequency Vector plot(f,fftshift(Pxx), 'linewidth', .25); axis([ -4 4 -80 0]) title('(c) Periodogram, Blackman ') ylabel('PSD, dB') set(gca, 'Xtick', -4: 2: 4) grid on xlabel('Frequency') subplot(224) NFFT = 256; x = data(1:NFFT); w = kaiser(length(x(1:NFFT)), 4) x = w.* x'; xf = fft(x,NFFT); Pxx=10*log10(xf.*conj(xf)); Pxx = Pxx - max(Pxx); f = (-NFFT/2)*fs/NFFT: fs/NFFT: (NFFT/2-1)*fs/NFFT; %Frequency Vector plot(f,fftshift(Pxx), 'linewidth', .25); axis([ -4 4 -80 0]) title('(d) Periodogram, Kaiser ') ylabel('PSD, dB') set(gca, 'Xtick', -4: 2: 4) grid on xlabel('Frequency') %% % Example from Poularkis book clf n = 1: 48; n2 = 1: 80; x = sin(.2*pi*n)+sin(0.22*pi*n)+ zeros(1,48); figure(1) subplot(211) grid on plot(n, x) title('Two close frequency sinusoids') ylabel('Amplitude') xlabel('Time') grid on subplot(212) px = 20*log10(fftshift(abs(fft(x, 512)))) px = px- max(px); plot(px) title('(Periodogram ') grid on axis([ 0 512 -40 0]) xlabel('Frequency') ylabel('PSD, dB') %% figure(2) y1 = [ x zeros(1, 32)] +.2*zeros(1,80); y2 = [zeros(1, 32) x] + .2*zeros(1,80); y = .2*y1 +((y2)*.2); px = 20*log10(abs(fft(y, 512))) px = px- max(px); figure(1) subplot(4,1,1) plot(n, x) title('(a) Signal with two sinusoids') ylabel('x') subplot(4,1,2) plot(n2, y1) title('(b) Part of signal with zer0-padding at end') ylabel('y1') subplot(4,1,3) plot(n2, y2) ylabel('y2') title('(c) Second part of signal with zer0-padding at start') subplot(4,1,4) plot(n2, y) ylabel('y1+y2') title('(d) Sum of both parts') xlabel('Time') figure(3) w = 0: 2*pi/512: 2*pi-2*pi/512; plot(w, px) ylabel('PSD') axis([ 0 6.3 -50 0]) grid on ylabel('PSD, dB') title('Periodogram of summed signal') xlabel('Frequency') % Chapter 7 %% clf N = 256; Fs = 1024; % fs time = 0: 1/Fs : 1; % one second of data x = 4*sin(2*pi*100*time) - 4/2*sin(2*pi*200*time) + 4/3*sin(2*pi*300*time) + randn(size(time)); pxx1 = fftshift(abs(fft(x(1:256), N).^2/256))... + fftshift(abs(fft(x(257:512), N).^2/256)) ... + fftshift(abs(fft(x(513:768), N).^2/256))... + fftshift(abs(fft(x(769:1024), N).^2/256)); pxx2 = fftshift(abs(fft(x(1:256), N).^2/256))... + fftshift(abs(fft(x(129:384), N).^2/256))... + fftshift(abs(fft(x(257:512), N).^2/256))... + fftshift(abs(fft(x(385:640), N).^2/256))... + fftshift(abs(fft(x(513:768), N).^2/256))... + fftshift(abs(fft(x(769:1024),N).^2/256))... + fftshift(abs(fft(x(641:896), N).^2/256)); wx = hanning(256)'; pxx3 = fftshift(abs(fft(wx.*x(1:256), N).^2/256))... + fftshift(abs(fft(wx.*x(129:384), N).^2/256))... + fftshift(abs(fft(wx.*x(257:512), N).^2/256))... + fftshift(abs(fft(wx.*x(385:640), N).^2/256))... + fftshift(abs(fft(wx.*x(513:768), N).^2/256))... + fftshift(abs(fft(wx.*x(641:896), N).^2/256))... + fftshift(abs(fft(wx.*x(769:1024), N).^2/256)); subplot(3,1,1) plot(((0:255)/256*Fs)- Fs/2, 10*log10(pxx1)) axis([ -500 500 -20 40]) title('(a) Averaged Periodogram, 4 Non-overlapping sgments') ylabel('PSD, dB') set(gca, 'Xtick', -500: 125: 500) grid on subplot(3,1,2) plot(((0:255)/256*Fs)- Fs/2, 10*log10(pxx2)) axis([ -500 500 -20 40]) title('(b) Averaged Periodogram, 7 Overlapping sgments') ylabel('PSD, dB') set(gca, 'Xtick', -500: 125: 500) grid on subplot(3,1,3) plot(((0:255)/256*Fs)- Fs/2, 10*log10(pxx3)) axis([ -500 500 -20 40]) grid on title('(c) Averaged Periodogram, 7 overlapping sgments with Hanning Window') grid on xlabel('Frequency') ylabel('PSD, dB') set(gca, 'Xtick', -500: 125: 500) % Chapter 8 %% % Welch Algorithm clf x = xlsread('nonlin',1, 'a1:a8192'); for numseg = 1:6; seglen = datalen/numseg; overf = .5; N = seglen; overlap = overf*seglen; NFFT = seglen; winx = hanning(seglen); seg1 = winx*x(1:seglen); X1 = (abs(fft(seg1,rectwin(seglen), N))); Xa = 10*log10((AllX )); maxA = max(Xa); n1 = 1: N/2+1; str = sprintf('Averged periodogram, no. of segment = %f', numseg); title(str); end seg1 = x(1:seglen); seg2 = x(overlap*1+1:overlap*3); seg3 = x(overlap*2+1:overlap*4) ; seg4 = x(overlap*3+1:overlap*5); seg5 = x(overlap*4+1:overlap*6) ; seg6 = x(overlap*5+1:overlap*7); seg7 = x(overlap*6+1:overlap*8); seg8 = x(overlap*7+1:overlap*9) ; X1 = (periodogram(seg1,rectwin(seglen), N)); X2 = (periodogram(seg2,rectwin(seglen), N)); X3 = (periodogram(seg3,rectwin(seglen), N)); X4 = (periodogram(seg4,rectwin(seglen), N)); X5 = (periodogram(seg5,rectwin(seglen), N)); X6 = (periodogram(seg6,rectwin(seglen), N)); X7 = (periodogram(seg7,rectwin(seglen), N)); X8 = (periodogram(seg8,rectwin(seglen), N)); AllX = (X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8)/8; Xa = 10*log10((AllX )); maxA = max(Xa); n1 = 1: N/2+1; str = sprintf('Averged periodogram, no. of segment = %f', numseg); title(str); xlabel('Frequency') ylabel('PSD, dB') n2 = -N/2: 0; figure(1) plot(n1, Xa - maxA, 'linewidth', .4) grid on axis([ -256 256 -30 0]) set(gca, 'Xtick', -datalen/2: datalen/8: datalen/2) hold on plot(fliplr(n2), Xa - maxA, 'linewidth', .4) hold on n3 = -N/2: N/2-1; plot((n3/2)-1, Xa) axis([ -datalen/2 datalen/2 -60 0]) %% % Chapter 8, Fig. 8.26 Convolution and Autocorrelation x1 = [ 1 1 -1 2 2 -1 -1 -1 1 2 2 -1 2] x2 = [ 1 1 2 2 -1 -1 -1 -1 -1 2 2 1 1 ] subplot(421) stem(x1) title('(a) A unsymmetrical sequence, s1') xlabel('Sample, n') ylabel('Amplitude') subplot(422) stem(x2) title('(a) A symmetrical sequence, s2') xlabel('Sample, n') subplot(425) axf1 = xcorr(x1, 'biased') axf1= axf1/max(axf1); stem(axf1) title('(c) ACF of s1') ylabel('ACF') subplot(426) axf2 = xcorr(x2, 'biased') axf2= axf2/max(axf2); stem(axf2) title('(c) ACF of s2') subplot(423) axf3 = conv(x1, x1); axf3= axf3/max(axf3); stem(axf3) title('(c) Self-convolution, s1') ylabel('Convolution') subplot(424) axf4 = conv(x2, x2); axf4= axf4/max(axf4); stem(axf4) title('(c) Self-convolution, s2') subplot(427) axf5 = conv(x1, flip(x1)) axf5= axf5/max(axf5); stem(axf5) title('(c) Self-convolution, s1, flipped s1') xlabel('Lag') ylabel('ACF') subplot(428) axf6 = conv(x2, x2) axf6= axf6/max(axf6); stem(axf6) title('(c) Self-convolution, s2, flipped s2') xlabel('Lag %% % ACF x = [ 1 2 1 4 3] [c, lags] = xcorr(x); y = xcorr(x); max(y) length(y) subplot(4,1,1) t = linspace( 1, length(y), 9); stem(lags, y, 'filled') ylabel('Raw') axis([-4 4 0 32]) y2 = xcorr(x, 'biased'); subplot(4,1,2) stem(lags, y2, 'filled'); axis([-4 4 0 10]) ylabel('Biased') y3 = xcorr(x, 'unbiased'); subplot(4,1,3) stem(lags, y3, 'filled'); axis([-4 4 0 10]) ylabel('Unbiased') max(y3) y4 = xcorr(x, 'coeff'); subplot(4,1,4) stem(lags, y4, 'filled'); axis([-4 4 0 2]) ylabel('Norm') xlabel('Lags') max(y3) %% x = [ 1 2 1.3 4 .3 .6]; x1 = [ 0 1 2 1.3 4 .3 .6]; x2 = [ 0 0 0 0 1 2 1.3 4 .3 .6]; x3 = [ 0 0 0 0 0 1 2 1.3 4 .3 .6]; n = 1: length(x3); subplot(4,1,1) stem(n(1:length(x)), x, 'filled') ylabel('x') axis([0 12 -1 5]) subplot(4,1,2) stem(n(2:7), x, 'filled') ylabel('lag = 1') axis([0 12 -1 5]) subplot(4,1,3) stem(n(4:9), x, 'filled') ylabel('lag = 4') axis([0 12 -1 5]) subplot(4,1,4) stem(n(6:11), x, 'filled') ylabel('lag = 6') axis([0 12 -1 5]) % Chapter 9 % Several periodograms % Welch clf x = xlsread('IQData',1, 'a1:a8192'); datalen = 8196; numseg = 8; overf = 1; seglen = datalen/numseg; N = seglen; overlap = overf*seglen; subplot(311) datalen = 8192; numseg = 1; overf = 1; seglen = datalen/numseg; N = seglen; overlap = overf*seglen; winx = hanning(seglen); seg1 = x(1:seglen).*winx; X1 = fftshift(abs(fft(seg(i), datalen))) Xa = 20*log10((X1 )); maxA = max(Xa); n1 = 1: N; plot(n1, Xa - maxA, 'linewidth', .4) axis([ 0 datalen -60 0]) ylabel('PSD, dB') str = sprintf('Number of segements %f',numseg); title(str); subplot(311) datalen = 8192; numseg = 4; overf = .5; seglen = datalen/numseg; N = seglen; overlap = overf*seglen; winx = hanning(seglen); seg2 = x(overlap*1+1:overlap*2); seg2 = x(overlap*1+1:overlap*3); seg3 = x(overlap*2+1:overlap*4) ; seg4 = x(overlap*3+1:overlap*5); seg1 = x(1:seglen).*winx; X1 = fftshift(abs(fft(seg1, datalen))) Xa = 20*log10((X1 )); maxA = max(Xa); n1 = 1: N; plot(n1, Xa - maxA, 'linewidth', .4) axis([ 0 datalen -60 0]) ylabel('PSD, dB') str = sprintf('Number of segements %f',numseg); title(str); %% clf x = xlsread('nonlin',1, 'a1:a8192'); datalen = 8192; N = 8192; n = -8192/2: 1: 8192/2 - 1 winx = hanning(seglen); numseg = 1 seglen = datalen/numseg; seg1 = x(1:seglen); xf =20*log10( fftshift(abs(fft(seg1, N)))); xfmax = max(xf); xf = xf - xfmax; %plot(n, xf) title('Segments = 1') axis([-8192/2 8192/2 -60 0]) subplot(221) numseg = 2 seglen = datalen/numseg; overlap = overf*seglen; seg1 = x(1:seglen) seg2 = x(1+overlap: seglen+overlap) xf =20*log10(fftshift(abs(fft(seg1, N)))+fftshift(abs(fft(seg2, N)))) xfmax = max(xf); xf = xf - xfmax; plot(n, xf) title('(a) Segments = 2') axis([-8192/2 8192/2 -60 0]) grid on subplot(222) numseg = 4 seglen = datalen/numseg overlap = overf*seglen; seg1 = x(1:seglen) seg2 = x(1*overlap+1: seglen+overlap) seg3 = x(2*overlap+1: seglen+2*overlap) seg4 = x(3*overlap+1: seglen+3*overlap) xf1 = fftshift(abs(fft(seg1, N))); xf2 = fftshift(abs(fft(seg2, N))); xf3 = fftshift(abs(fft(seg3, N))); xf4 = fftshift(abs(fft(seg4, N))); allxf = xf1 + xf2 + xf3 + xf4; allxf = 20*log10(allxf) xfmax = max(allxf); xf = allxf - xfmax; plot(n, xf) title('(b) Segments = 4') axis([-8192/2 8192/2 -60 0]) grid on subplot(223) numseg = 8 seglen = datalen/numseg overlap = overf*seglen; seg1 = x(1:seglen) seg2 = x(1*overlap +1: seglen+overlap) seg3 = x(2*overlap +1: seglen+2*overlap) seg4 = x(3*overlap +1: seglen+3*overlap) seg5 = x(4*overlap +1: seglen+4*overlap) seg6 = x(5*overlap +1: seglen+5*overlap) seg7 = x(6*overlap +1: seglen+6*overlap) seg8 = x(6*overlap +1: seglen+7*overlap) xf1 = fftshift(abs(fft(seg1, N))); xf2 = fftshift(abs(fft(seg2, N))); xf3 = fftshift(abs(fft(seg3, N))); xf4 = fftshift(abs(fft(seg4, N))); xf5 = fftshift(abs(fft(seg5, N))); xf6 = fftshift(abs(fft(seg6, N))); xf7 = fftshift(abs(fft(seg7, N))); xf8 = fftshift(abs(fft(seg8, N))); allxf = xf1 + xf2 + xf3 + xf4 + xf5 + xf6 + xf7 + xf8; allxf = 20*log10(allxf) xfmax = max(allxf); xf = allxf - xfmax; plot(n, xf) title('(c) Segments = 8') axis([-8192/2 8192/2 -60 0]) grid on subplot(224) numseg = 16 seglen = datalen/numseg overlap = overf*seglen; seg1 = x(1:seglen) seg2 = x(1*overlap +1: seglen+overlap) seg3 = x(2*overlap +1: seglen+2*overlap) seg4 = x(3*overlap +1: seglen+3*overlap) seg5 = x(4*overlap +1: seglen+4*overlap) seg6 = x(5*overlap +1: seglen+5*overlap) seg7 = x(6*overlap +1: seglen+6*overlap) seg8 = x(7*overlap +1: seglen+7*overlap) seg9 = x(8*overlap +1: seglen+8*overlap) seg10 = x(9*overlap +1: seglen+9*overlap) seg11 = x(10*overlap +1: seglen+10*overlap) seg12 = x(11*overlap +1: seglen+11*overlap) seg13 = x(12*overlap +1: seglen+12*overlap) seg14 = x(13*overlap +1: seglen+13*overlap) seg15 = x(14*overlap +1: seglen+14*overlap) seg16 = x(15*overlap +1: seglen+15*overlap) xf1 = fftshift(abs(fft(seg1, N))); xf2 = fftshift(abs(fft(seg2, N))); xf3 = fftshift(abs(fft(seg3, N))); xf4 = fftshift(abs(fft(seg4, N))); xf5 = fftshift(abs(fft(seg5, N))); xf6 = fftshift(abs(fft(seg6, N))); xf7 = fftshift(abs(fft(seg7, N))); xf8 = fftshift(abs(fft(seg8, N))); xf9 = fftshift(abs(fft(seg9, N))); xf10 = fftshift(abs(fft(seg10, N))); xf11 = fftshift(abs(fft(seg11, N))); xf12 = fftshift(abs(fft(seg12, N))); xf13 = fftshift(abs(fft(seg13, N))); xf14 = fftshift(abs(fft(seg14, N))); xf15 = fftshift(abs(fft(seg15, N))); xf16 = fftshift(abs(fft(seg16, N))); allxf = xf1 + xf2 + xf3 + xf4 + xf5 + xf6 + xf7 + xf8 + xf9 + xf10 + xf11 + xf12 + xf13 + xf14 + xf15 + xf16; allxf = 20*log10(allxf) xfmax = max(allxf); xf = allxf - xfmax; plot(n, xf) title('(d) Segments = 16') axis([-8192/2 8192/2 -60 0]) grid on
2 comments on “Matlab code for FFT book
  1. Sven E Widmalm says:

    You wrote: Posted on November 7, 2016 by Charan L. — No Comments ↓
    Intuitive Guide to Fourier Analysis and Spectral Estimation
    Matlab code for selected figures and examples is given here.
    If you find any issues with these programs, please email me.
    Last edited 11/24/2017
    Matlab files from the book are included here.
    Signal data required for some of these scripts will be uploaded soon.

    What is your definition of “soon”?
    Sven

Leave a Reply

Your email address will not be published.

*

This site uses Akismet to reduce spam. Learn how your comment data is processed.