clear
%PSD Test Program Demonstrating Concepts with a White Noise Spectrum
%set a sampling interval (s)
delta=1;

%sampling rate (Hz)
fs=1/delta;

%data length in samples
N=2^14;
t = (0:N-1)*delta;

%synthetic white noise signal, units U
x=(1e-6)*randn(N,1);
figure(1)
plot(t,x)
title ('Time Series')
ylabel('U')
xlabel('t (s)')

%synthetic white noise signal with a 1 Hz sinusoid added for demo purposes
%x=(1e-6)*(randn(N,1)+sin(2*pi*1*(0:N-1)'*delta));

disp('Time Domain Estimates')

%integrated signal energy
p1=sum(x.*x)*delta;
disp(['total signal energy: ',num2str(p1),' (U^2)*s'])

%signal power
p2=p1/(N*delta);
disp(['signal power:  ',num2str(p2),' (U^2)'])

%power spectral density for a white process is equally distributed over the Nyquist interval
p3=p2/(fs/2);
disp(['1-sided signal power per Hz:  ',num2str(p3),' (U^2/Hz)'])
p3db=10*log10(p3);
disp(['1-sided signal power per Hz:  ',num2str(p3db),' (dB rel. 1 U^2/Hz)'])

disp(sprintf('\n'))
disp('Frequency Domain Estimates')

[Pxx,F]=pmtm(x,4,[],fs);
mPxx=mean(Pxx);
Pxxdb=10*log10(mPxx);
figure(2)
semilogx(F,10*log10(Pxx),'g');
hold on
semilogx(F,Pxxdb*ones(N/2+1,1),'r');
hold off
title ('1-sided time series PSD')
ylabel ('dB rel. 1 U^2/Hz')
xlabel ('Hz')

p4=sum(Pxx)*(F(2)-F(1));

%total energy 
disp(['signal energy estimated from PSD:  ',num2str(p4*N*delta),' (U^2*s)'])

%total power integrated from the PSD estimate
disp(['integrated PSD (power):  ',num2str(p4),' (U^2)'])

%power spectral density mean level
disp(['PSD mean PSD level (power/Hz):  ',num2str(mPxx),' (U^2/Hz)'])
disp(['PSD mean 1-sided signal power per Hz:  ',num2str(Pxxdb),' (dB rel. 1 U^2/Hz)'])