代码如下;
%u是导入的参数3600*1的矩阵
Fs=1;
nfft=2^nextpow2(length(u));
L = size(u,1); % 信号长度
F = Fs/2*linspace(0,1,nfft/2+1); %(0:1:N/2)*Fs/N
Y = fft(u,nfft)/L; %除以L乘以2
y = 2*abs(Y(1:nfft/2+1,:));
window=hanning(nfft)
noverlap=nfft/2;
[Sv,f]=pwelch(u,window,noverlap,nfft,Fs,'onesided')
fn=z*f/v0;
Sv=Sv.*f/sdz;