随机波浪的谱分析

随机波浪的谱分析

常用的二维波能谱分析方法有两种，一种为自相关函数（协方差函数）法，另一种为快速傅里叶变换法。自相关函数法含有能量分布的信息，快速傅里叶变换法则更加直接快速。

自相关函数法

自相关函数法主要思路是首先求解自相关函数，然后对自相关函数进行傅里叶变换得出波浪谱的粗谱，最后通过窗函数对数据进行平滑。

1. 求自相关函数

假设现有M个规则波叠加产生的随机波浪可表示为：
$$\zeta (x,t)=\sum _{i=1}^M{a}_i\cos (k_{i}x-{\omega }_{i}t+{\epsilon }_i)$$
$\zeta (x,t)$为齐次平稳的，具有各态历经性，数学期望$\mathrm{E}[\zeta (x,t)]=0$，自相关函数如下式表示：
$$R_{\eta \eta }(\tau )=\mathrm{E}[\zeta (x,t)\zeta (x,t+\tau )]=\overline{\zeta (x,t)\zeta (x,t+\tau )}$$
上面这个是数学公式，但是在实际操作中由于采样数据为离散的，所以要对上式进行离散化处理，方法如下：
假设采样数据点数为N，采样时间间隔为$\Delta t$，则相关函数可以表示为
$$\begin{gathered} R(\tau )=R(\nu \Delta t)=\frac{1}{N-\nu }\sum _{n=1}^{N-\nu }\zeta (t_n+\nu \Delta t)\zeta (t_n) \\ \tau =\nu \Delta t,\nu =0,1,2\dots m \end{gathered}$$
这其中的m与前文中所述M是可以理解为同一个值，这是因为靶谱中离散设定的频率个数是M个，那么理论上讲频谱中也就应该有M个。当然我们在开始时并不知道造波时采用了多少个离散频率在波浪中，所以可以根据情况进行设定，在后面做谱的时候，设定的m为多大，就会有多少个离散频率值，设的m值越大，谱中的频率数越多。

将上式拆开进行分析，如下：
$$\begin{gathered} \nu =0R(0)=\frac{1}{N}[\zeta (t_1)\zeta (t_1)+\zeta (t_2)\zeta (t_2)+\dots +\zeta (t_{N-1})\zeta (t_{N-1})+\zeta (t_N)\zeta (t_N)] \\ \nu =1R(\Delta t)=\frac{1}{N-1}[\zeta (t_1+\Delta t)\zeta (t_1)+\zeta (t_2+\Delta t)\zeta (t_2)+\dots +\zeta (t_{N-2}+\Delta t)\zeta (t_{N-2})+\zeta (t_{N-1}+\Delta t)\zeta (t_{N-1})] \\ \nu =2R(2\Delta t)=\frac{1}{N-2}[\zeta (t_1+2\Delta t)\zeta (t_1)+\zeta (t_2+2\Delta t)\zeta (t_2)+\dots +\zeta (t_{N-3}+2\Delta t)\zeta (t_{N-3})+\zeta (t_{N-2}+2\Delta t)\zeta (t_{N-2})] \\ \dots \\ \nu =mR(m\Delta t)=\frac{1}{N-m}[\zeta (t_1+m\Delta t)\zeta (t_1)+\zeta (t_2+m\Delta t)\zeta (t_2)+\dots +\zeta (t_{N-m}+m\Delta t)\zeta (t_{N-m})] \end{gathered}$$

上式中每一个$R(\nu \Delta t)$均可以求出相应的值，编程中仅需采用一次循环就可实现，matlab代码如下：

Data=Surf-mean(Surf);%变量Surf为波面采样数组，变量Data是变量Surf数学期望归零的数组
N=length(Data);%获得采样点个数
dt=1/fs;%此处fs为采样频率，可以在前期前面进行设定，dt为采样周期
M=N/30;%这里设置M为采样点个数的1/30个

% 自相关方程计算
for v=0:M
    R(v+1)=sum(Data(1+v:N).*Data(1:N-v))/(N-v);
end

2. 进行傅里叶变换，求得波浪谱粗谱

自相关函数与谱密度函数存在Fourier变换关系，即所谓的纳维-辛钦定理，公式如下：
$$\begin{gathered} S(\omega )=\frac{2}{\pi }{\int }_0^{\infty }R(\tau )\cos (\omega \tau )d\tau \\ 或 \\ S(f)=4{\int }_0^{\infty }R(\tau )\cos (2\pi f\tau )d\tau \end{gathered}$$
谱密度函数分布的频率范围为$[0,f_N]$，其中$f_N$为奈奎斯特频率（Nyquist frequency)，等于离散信号采样频率的一半，为
$$f_N=\frac{{\omega }_N}{2\pi }=\frac{1}{2\Delta t}$$
令$f_n=\frac{n}{m}{f}_N,n=0,1,2\dots m$，则
$$S^{\prime }(f_n)=\frac{2}{\pi }\sum _{\tau =0}^{m\Delta t}R(\tau )\cos ({\omega }_n\tau )d\tau =\frac{2}{\pi }\sum _{\nu =0}^{m}R(\nu \Delta t)\cos ({\omega }_n\nu \Delta t)\Delta t$$
考虑到${\omega }_n=2\pi f_n=2\pi \frac{n}{m}{f}_N=2\pi \frac{n}{m}\frac{1}{2\Delta t}$，则${\omega }_n\Delta t=\frac{\pi n}{m}$,并采用梯形法求积分，得到如下式子：
$$\begin{gathered} S^{\prime }(f_n)=\frac{2\Delta t}{\pi }[\frac{1}{2}R(0)+\sum _{\nu =1}^{m-1}R(\nu \Delta t)\cos \frac{\pi \nu n}{m}+\frac{1}{2}R(m\Delta t)\cos (n\pi )] \\ n=0,1,2,\dots ,m \end{gathered}$$
matlab代码实现如下：

df=1/(2*M*dt);
for n=0:M
    L_sum=0;
    for v=1:M-1
        L_sum=R(v+1)*cos(2*pi*(n-1)*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

3. 对粗谱进行光顺

以上代码中的L(N)即为粗谱，但该谱显示出来可能会如下图所示：

所以，此时需要对粗谱进行光顺。一般采用光顺的方法是采用窗函数，所谓的窗函数就是采用特定的系数对谱中的数据前后相邻数值进行加权平均，第一种采用的是hamming window, 如下所示：
$$S(f_n)=\{\begin{array}{l}0.54S^{\prime }(f_0)+0.46S^{\prime }(f_1)n=0\\ 0.23S^{\prime }(f_{n-1})+0.54S^{\prime }(f_n)+0.23S^{\prime }(f_{n+1})n=1,2,\dots ,m-1\\ 0.46S^{\prime }(f_{n-1})+0.54S^{\prime }(f_n)n=m\end{array}$$
第二种为hanning window法：
$$S(f_n)=\{\begin{array}{l}0.5S^{\prime }(f_0)+0.5S^{\prime }(f_1)n=0\\ 0.25S^{\prime }(f_{n-1})+0.5S^{\prime }(f_n)+0.25S^{\prime }(f_{n+1})n=1,2,\dots ,m-1\\ 0.5S^{\prime }(f_{n-1})+0.5S^{\prime }(f_n)n=m\end{array}$$
matlab实现代码如下：

%hamming window function
    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 window function
    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;
    Frequency=x*2*pi;

参考资料：
随机波浪谱分析——自相关函数法 （MATLAB）.
浮子式波浪能转换装置机理的频域及时域研究

快速傅里叶变换法

快速傅里叶变换法是通过对波浪时历曲线进行快速傅里叶变换，然后再将数据转换成谱数据得到谱图。以往关于快速傅里叶变换法得到的波浪谱的公式推导都未有比较明确的描述，导致我们并不清楚最后的公式如何运用，此处我将对此方法进行公式推导，以减少读者困惑。
$$\begin{gathered} a_i=\sqrt{2S_{\eta \eta }({\hat{\omega }}_i)\Delta {\omega }_i}⟹S_{\eta \eta }({\hat{\omega }}_i)=\frac{a_i^2}{2\Delta \omega } \\ \because \Delta \omega =2\pi \Delta f,\Delta f=\frac{Fs}{N},Fs=\frac{1}{\Delta T} \\ \therefore \Delta \omega =\frac{2\pi }{N\cdot \Delta T} \\ ⟹S_{\eta \eta }({\hat{\omega }}_i)=\frac{a_i^{2}N\cdot \Delta T}{4\pi } \end{gathered}$$
$W(\Delta T\cdot i)$为波浪采样点，$N$为采样点数量，$\Delta T$为采样周期，
$$\begin{gathered} Y_i=abs[FFT(W(\Delta T\cdot i))],则 \\ {a}_i=\frac{2Y_i}{N}，当i=2\sim N/2时; \\ \therefore S_{\eta \eta }({\hat{\omega }}_i)=\frac{\Delta T}{\pi \cdot N}\cdot Y_i^2,当i=2\sim N/2时. \end{gathered}$$
以上这种形式就和在文献中看到的相同了。
由于在matlab快速傅里叶变换中$Y_1$是零频时（即波浪的直流分量）的值，而波浪是经过归零处理的，即零频时的值应为0，所以$a_1$的值应为零，与上式统一起来就是
$$S_{\eta \eta }({\hat{\omega }}_i)=\frac{\Delta T}{\pi \cdot N}\cdot Y_i^2，i=1\cdots \frac{N}{2}$$
具体文献可参看快速傅里叶变换法的海浪谱估计

matlab实现方式

%% Spectrum Analysis
%   Code by YY-ZHU
%   Modified on 6/3/2021
%% 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,Frequency]=Spectrum_Analysis_FFT(b,Fs,Method)
  T = 1/Fs;             % Sampling period       
L = length(b);             % Length of signal
if mod(L,2)~=0
    a=b(1:end-1);
    L=L-1;
else
    a=b;
end

Y = fft(a);
P2 = abs(Y/L);  % 2-sides spectrum
P1 = P2(1:L/2+1); % fetch half of P2
P1(2:end-1) = 2*P1(2:end-1); 

for i=1:L/2+1
    S2(i)=(P1(i)^2)*T*L/4/pi;
    %S2(i)=1/(2*pi*L)*P1(i)^2;
end
for i=1:100
S2=WindowFunction(S2,Method);
end
Frequency = 2*pi*Fs*(0:(L/2))/L;
Spectrum = S2;

其中window function如下

%% Information
%   this function include two methods to smooth the spectrum
%   1. Hammming window
%   2. Hanning window
%   the second argument is set as 'Hammming' or 'Hanning' to select the
%   method

function [S1]=WindowFunction(P1,Method)
L=length(P1);

if  Method == ['Hamming']
    for i=1:L-2
        S1(i+1)=0.23*P1(i)+0.54*P1(i+1)+0.23*P1(i+2);
    end
    S1(1)=0.54*P1(1)+0.46*P1(2);
    S1(end)=0.46*P1(end-1)+0.54*P1(end);
elseif Method == ['Hanning']
    for i=1:L-2
        S1(i+1)=0.25*P1(i)+0.5*P1(i+1)+0.25*P1(i+2);
    end
    S1(1)=0.5*P1(1)+0.5*P1(2);
    S1(end)=0.5*P1(end-1)+0.5*P1(end);
else
    disp('wrong input');
    return;
end