现在测得波面位移随时间的变化，有很多数据，下面的程序可以给出海浪谱。

wave=wave2(3000:end); %波面位移测量从第3000个开始，根据自己的实际情况修改
L=length(wave);
Fs=50;  %采样频率，根据自己的实际情况修改
T=1/50; %根据自己的实际情况修改
NFFT = 2^(nextpow2(L)-1);
y = fft(wave,NFFT)/(NFFT/2);
f = Fs*(0:NFFT/2)/NFFT; %只画频率小于折叠频率的部分
zfu=abs(y(1:NFFT/2+1)); %振幅
Sw=(zfu.^2)/2;  %能量
sigma1=sum(Sw(2:end)) %方差，Sw(1)异常大，对应的是频率为0的振幅，不参加后面的运算
%plot(f(2:end),Sw(2:end)/(f(2)-f(1)))
ff=(2*pi*f).^2;
Sw2=ff'.*Sw;
sigma2=sum(Sw2(2:end)); %谱的二阶矩
H_avg=sqrt(2*pi*sigma1)                 %求得平均波高为0.9801m
T_avg=2*pi*sqrt(sigma1/sigma2)      %求得平均周期为0.2890s
%画出平滑处理后的海浪谱
plot(f(2:end),smooth(Sw(2:end)/(f(2)-f(1)),50,'lowess'))
xlabel('w')
ylabel('S(w）')

http://s13/mw690/003cvAxDzy6WKeb4CcIfc&690

使用上跨零法统计分析程序如下：

count=1;
num=size(wave);
for i=1:num-1
    if wave(i)<0 &&wave(i+1)>0
        up(count)=i;
        count=count+1;
    end
end
num=size(up)
for j=1:num(2)-1
    hh=max(wave(up(j):up(j+1)));
    ll=min(wave(up(j):up(j+1)));
    H(j)=hh-ll;
    T(j)=(up(j+1)-up(j))/50;
end
%采用上跨零法
mean(H) %求得平均波高1.0588m
mean(T) %求得平均周期0.3220s

我们谱分析得到的平均波高略小于统计分析得到的，这是自然的，因为只算了折叠频率以下的部分；0.322/0.289=1.1142,小于理论值1.2，原因不明，可能是误差所致。