For a project I’m currently working on, I composed a irregular ocean wave (function of time) by the summation of n regular waves. For the validation and further processing I’d like to produce the following graphs:
- Power Spectral Density (PSD) spectrum of the wave (single sided) with:
- Units: x-axis: [Hz], y-axis: [m^2/Hz]
- Units: x-axis: [rad/s], y-axis: [m^2/(rad/sec)]
Question 2 In order to produce the single sided PSD graphs I found some code online. The peaks are shown at the right period or frequency however, the information on the y-axis has to be scaled some way (changing the amount of samples per second (i.e. parameter "Fs") or lenght of the signal (i.e. parameter "StopTime") results in other values on the y-axis). Is this value shown in "bins"? How to change this correctly when the periodogram function is used?
CODE 2
Fs = 6; % samples per second
dt = 1/Fs; % seconds per sample
StopTime = 30*3600; % seconds
t = (0:dt:StopTime-dt)'; % seconds
N = size(t,1);
W = 5*cos(t) + 10*cos(8*t);
p = 5; %number of periodograms
%p=10;
%p=5;
%W is divided into 20 samples/periodograms 18000/20=900 samples/periodogram
%W is divided into 10 periodograms 18000/10=1800 samples/peridogram
%W is divided into 5 periodograms 18000/5=3600 samples/periodogram
%Very long records needed!!! 90 min for 10 periodograms or 45 for 5.
%Periodogram samples
%x=[];
for i=1:p
x(:,i)=W((1+(i-1)*N/p):((N/p)*i));
end
%Next power of 2 greater or equal to length of the sample vector to
%calculate fft. FFT needs a length power of 2 to be efficient. For a
%shorter record, zeros will be added to reach this value.
nfft=2^(nextpow2(length(x(:,1))));
%Calculate PSD
Pxx = [];
for i=1:p
[Pxx(:,i),f]=periodogram(x(:,i),[],nfft,Fs);
end
mW=sum(Pxx,2)/p; %Average periodograms
ft=1./f; %Wave period vector for the plot. f represents the frequency vector
%Wave frequency spectral density plot
figure;
plot(f,mW,'-')
xlabel('Frequency (Hz)')
title('Power wave Spectrum')
grid;
%Wave Spectral density plot with wave period
figure;
plot(ft,mW);
xlabel('Wave period (s)')
xlim([0 25])
title('Power wave Spectrum')
grid;