0
$\begingroup$

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;
$\endgroup$

migrated from stackoverflow.com Nov 20 '14 at 13:54

This question came from our site for professional and enthusiast programmers.

  • 1
    $\begingroup$ in which way are your amplitudes incorrect? Please describe what you get and what you expect instead. $\endgroup$ – SBH Nov 19 '14 at 15:17
  • $\begingroup$ @SBH: Regarding to Question 1 and CODE 1: The amplitudes are respectively 4.4 and 9.72 (determined reading the graph). The input signal has amplitudes of 5 and 10, the time signal is long enough and the stepsize small enough. I expect values of 5 and 10. Regarding to question 2: In my opinion it's incorrect that the PSD changes when the lenght of the signal or step size changes. $\endgroup$ – knarf Nov 19 '14 at 16:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy