随机波浪谱分析——自相关函数法 （MATLAB）

基本原理

自相关函数

假设采样时距为$\Delta t$，则海浪随机序列自相关函数可以写成
$R(\upsilon \Delta t)=\frac{1}{N-\upsilon }\sum _{n=1}^{N-\upsilon }x(t_n+\upsilon \Delta t)x(t_n)$
$\tau =\upsilon \Delta t,\upsilon =0,1,2,...,m$
这样可以得到$R(\upsilon )$的m+1个值，他们具有等间隔$\Delta t$。

粗谱

这里我们令$L_n$代表频率$f_n$对应的谱初值，有
$L_n=\frac{2}{\pi }\sum _{\tau =0}^{m\Delta t}R(\tau )\cos {\omega }_{\mathrm{n}}\tau \mathrm{d}\tau$

谱平滑

谱的平滑可采用哈明窗（Hamming）或哈宁窗（Hanning）函数进行改进

哈明窗

$S({\omega }_n)=0.23L_{n-1}+0.54L_n+0.23L_{n+1}$
对于两个端点
$S({\omega }_0)=0.54L_0+0.46L_1$
$S({\omega }_n)=0.46L_{m-1}+0.54L_m$

哈宁窗

$S({\omega }_n)=0.25L_{n-1}+0.5L_n+0.25L_{n+1}$
对于两个端点，系数均为0.5。

%% Spectrum Analysis
%   Code by JN_Cui
%   Modified on 23/6/2020
%% Defination
%   Three methods are included:
%   1. Self-Correlation Function (SCF) method
%   2. Fast Fourier Transform (FFT) method
%   3. Maximum Entropy Method (MEM)
%   The function for these three methods are named as:
%   1. Spectrum_Analysis_SCF
%   2. Spectrum_Analysis_FFT
%   3. Spectrum_Analysis_MEM
%% Function Spectrum_Analysis_SCF
function [Spectrum,Frequency]=Spectrum_Analysis_SCF(Surf,fs,Method)
Data=Surf-mean(Surf);
N=length(Data);
dt=1/fs;
M=roud(N/30);
df=1/(2*M*dt);
% Self-correlation function
for v=0:M
    R(v+1)=sum(Data(1+v:N).*Data(1:N-v))/(N-v);
end
% Spectrum
for n=0:M
    L_sum=0;
    for v=1:M-1
        L_sum=R(v+1)*cos(2*pi*(n)*df*v*dt)+L_sum;
    end
    L(n+1)=2*dt/pi*(R(1)/2+L_sum+1/2*R(M+1)*cos(2*pi*n*df*M*dt));
end
% Smooth the spectrum by Hamming
if Method == ['Hamming']
    for i=1:M-1
        S1(i+1)=0.23*L(i)+0.54*L(i+1)+0.23*L(i+2);
    end
    S1(1)=0.54*L(1)+0.46*L(2);
    S1(M+1)=0.46*L(M)+0.54*L(M+1);
    x=0:1/(2*M*dt):1/(2*dt);
    Spectrum=S1;
    Frequency=x*2*pi;
    % Hanning
elseif Method == ['Hanning']
    for i=1:M-1
        S2(i+1)=0.25*L(i)+0.5*L(i+1)+0.25*L(i+2);
    end
    S2(1)=0.5*L(1)+0.5*L(2);
    S2(M+1)=0.5*L(M)+0.5*L(M+1);
    x=0:1/(2*M*dt):1/(2*dt);
    Spectrum=S2(2:end);
    Frequency=x(1:end-1)*2*pi;
else
    disp('wrong input')
    Spectrum=NaN;
    Frequency=NaN;
    return
end