% 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).
%
%
clear
%Set Signal Parameters
fs = 8000;  %  Sampling frequency in Hz
sdur = 4;   % Signal duration in seconds
f0 = 1000;   %  Signal Frequency
snrdb = -10;  %  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)
%  Hamming is a tapering window on each data segment, try "boxcar(wl)" to see
%  effect from no tapering window
tapering = hamming(wl);  %  Create tapering window

%  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
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(['Comparing sinudoidal signal with and without additive WGN at SNR = ' num2str(snrdb)] )
legend({'With Noise', 'Original'})

%  Now compute PSD using the pwelch command

[p,f] = pwelch(r,tapering,wolap,nfft,fs);
figure
%  Plot in kHz
plot(f/1000,10*log10(p),'k')
xlabel('kHz')
ylabel('PSD')
title(['PSD of signal with noise at SNR = ' num2str(snrdb)])
%  Uncomment next set of lines to overlay signal PSD without noise
% hold on
% [p,f] = pwelch(s,tapering,wolap,nfft,fs);
% plot(f/1000,10*log10(p),'r')
% hold off
% title(['PSD of signal with and without noise at SNR = ' num2str(snrdb)] )
% legend({'With Noise', 'Original'})


