% 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
% Updated Sept. 2015
%
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'})