% Power Spectral Density Example using Pwelch.
% The script generates a sinusoid, normalizes its power, adds
% white Gaussian noise (WGN) at a specified SNR level, and
% then computes and plots the power spectral density (PSD).
%
%   Written by Kevin D. Donohue (donohue@engr.uky.edu) Sept. 2007
%
clear
%Set Signal Parameters
fs = 8000;  %  Sampling frequency in Hz
sdur = 4;   % Signal duration in seconds
f0 = 1000;   %  Signal Frequency
snrdb = 12;  %  Signal to noise ratio
%  PSD Parameters
wl = 512;  %  Data window length
wolap = fix(wl/2);  %  Number of overlapping points between segments
nfft = 2*wl;  %  Number of FFT points (either truncates or zero padds)


%  Create time axis
t = (0:round(fs*sdur)-1)/fs;
%  Create signal with unit power (RMS = 1) 
s = cos(2*pi*f0*t);
s = s/sqrt(mean(s.^2));  %  Normalize by RMS value (for sinusoid can also divide by Amplitude/sqrt(2))
%  Create zero mean Gaussian white noise with unit power
n = randn(size(s));  % Make noise vector same size as signal
%  Compute scaling factor for SNR
nlev = 10^(-snrdb/20);  %  Convert SNR from dB to linear scaling
%  Add noise to signal at specified SNR level
r = s+nlev*n;
%  Plot signals over a single window just to check things out
figure(1)
plot(t(1:wl),r(1:wl),'k',t(1:wl),s(1:wl),'r')  %  Plot signal with corrupted noise version
xlabel('Seconds')
ylabel('Amplitude')
title('Original signal (in red) Noise corrupted (in black)')

%  Now compute PSD using the pwelch command
%  Hamming is a tapering window on each data segment, try "boxcar(wl)" to see
%  effect from no tapering window
[p,f] = pwelch(r,hamming(wl),wolap,nfft,fs);
figure(2)
%  Plot in kHz, scale psd value by delta F to get power values
plot(f/1000,10*log10(p*fs/wl),'k')
xlabel('kHz')
ylabel('PSD')
title('PSD of of noise corrupted signal (in black)')
%  Uncomment next set of lines to the end to overlay signal PSD
%  without noise
% hold on
% [p,f] = pwelch(s,hamming(wl),wolap,nfft,fs);
% plot(f/1000,10*log10(p*fs/wl),'r')
% hold off
% title('PSD of Original signal (in red) Noise corrupted (in black)')