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