% Copyleft (C) Forrest Sheng Bao 1984-2009
% with Dept. of Electrical and Computer Engineering
% and  Dept. of Computer Science
% at Texas Tech University, Lubbock, Texas, USA
% http://fsbao.net
%
% This program compares the difference between Power Spectrum Density (PSD)
% of direct FFT and the PSD of autocorrelation function
%
% This program is a free software.
% It is licensed in GNU GENERAL PUBLIC LICENSE (GPL) v3 or later.
% So this program comes with ABSOLUTELY NO WARRANTY.
%
% You have the rights to freely COPY, MODIFY and REDISTRIBUTE this program
% under the terms of GNU GPL v3 or later, including KEEP the modified
% and/or redistributed code in GNU GPL v3 or later.
%
% If you do not know what GNU GPL is, please visit
% www.gnu.org/copyleft/gpl.html

clear all
close all

Fs = 1000;%sampling frequency

y = sig_noise([200], [0], 512);% generating the signal
L = length(y);%length of the signal

% Job 1: plot time domain signal
figure('Position',[100 100 900 700]);
subplot(2,2,1)
plot(1:L,y);
axis tight
axis 'auto y'
title('512-point signal containing one sine at 200 Hz, with SNR of 0dB')
xlabel('time (milliseconds)')
ylabel('amplitude')

% Job 2: do direct FFT
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);%only one side of the spectrum

% Plot single-sided amplitude spectrum.
% figure();
subplot(2,2,2)
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Single-sided amplitude spectrum of direct FFT')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

% Job 3: compute autocorrelation
y2 = xcorr(y,'coeff');

% plot autocorrelation sequence
% figure();
subplot(2,2,3)
plot(1:length(y2), y2)
axis tight
%axis 'auto y'
title({'Since the signal has high randomness,' 'it has only one autocorrelation peak when the shift equals to the signal length '})
xlabel('convolution shift (sample)')
ylabel('normalized autocorrelation value')

% Job 4: do FFT on autocorrelation sequence
L = length(y2);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y2,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);%only one side of the spectrum

% Plot single-sided amplitude spectrum.
%figure();
subplot(2,2,4)
plot(f,2*abs(Y(1:NFFT/2+1)))
title('PSD of the autocorrelation function')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

Published with MATLAB® 7.8