Fs = 6095;
x=double(squeeze(kh(301,403,15,:)));
% 各种谱估计算法
reply = input('Please Input spectrum estimation type:[1-9] ', 's');
switch str2double(reply)
case 1
Hs=spectrum.periodogram; %周期图法
figure;psd(Hs,x,'Fs',Fs)
case 2
Hs=spectrum.welch; %welch算法
figure;psd(Hs,x,'Fs',Fs)
case 3
Hs=spectrum.yulear;
figure;psd(Hs,x,'Fs',Fs)
case 4 %Thomson multitaper spectrum,汤姆森多窗口谱法
nw = 4; %控制分辨率，使用的斯莱皮恩蜡烛（Slepian tapers）宽度为2*nw-1
nfft = max(256,2^nextpow2(length(x)));
figure;pmtm(x,nw,nfft,Fs)
case 5
Hs=spectrum.mcov; %Modified covariance spectrum,修改后的协方差谱
figure;psd(Hs,x,'Fs',Fs)
case 6
Hs=spectrum.cov; %covariance spectrum,协方差谱
figure;psd(Hs,x,'Fs',Fs)
case 7
[a,e,k] = arburg(x,96); %求反射系数
a = abs(a);
order = numel(a(a>mean(a)))*2; %根据反射系数确定阶数选择
nfft = 256;
figure;pburg(x,order,nfft,Fs)
case 8 % MUSIC算法
segLen = 32;
ovlpPct = 95; %overlap百分比
Noverlap = ceil(ovlpPct/100*segLen); %overlap点数
nfft = max(256,2^nextpow2(segLen)); %FFT点数
thresh = 5; %噪声子空间门限值=λmin*thresh
nsinusoids = 6*2*2;%信号子空间维数,如果是实信号，维数要加倍。
P = [nsinusoids thresh];
win = hamming(segLen); %对输入数据进行分割加窗处理
figure;pmusic(x,P,nfft,Fs,win,Noverlap);
case 9 % Eigenvector spectrum 特征矢量谱算法（特征向量法）
segLen = 32;
ovlpPct = 95; %overlap百分比
Noverlap = ceil(ovlpPct/100*segLen); %overlap点数
nfft = max(256,2^nextpow2(segLen)); %FFT点数
thresh = 5; %噪声子空间门限值=λmin*thresh
nsinusoids = 6*2*4;%信号子空间维数,如果是实信号，维数要加倍。
P = [nsinusoids thresh];
win = hamming(segLen); %对输入数据进行分割加窗处理
figure;peig(x,P,nfft,Fs,win,Noverlap);
otherwise
warning(msgID,'Unexpected reply type. No plot created.');
end
%%
fs=1000;
u0av=zeros(17473);
u0av=u0av(:,1);
v0av=zeros(17473);
v0av=v0av(:,1);
for i=1:17473
u0av(i)=sum(u0(1:31,i))/31;
v0av(i)=sum(v0(1:31,i))/31;
end
xn1=u0av;
xn2=v0av;
%%
nfft=18000;
window=boxcar(length(xn1));
p1=floor(length(xn1)/3)+1;
[px1,f1]=pyulear(xn1,120,nfft,fs);%500<length(xn)
[px2,f2]=pyulear(xn2,120,nfft,fs);
px1max=max(px1);
px1=px1/px1max;
px2max=max(px2);
px2=px2/px2max;
for i=1:9001
x(i)=1000/f1(i);
end
x=log10(x);
figure('position',[0 0 700 215]);
set(gca,'position',[0.09,0.23,0.89,0.72]);
plot(x,px1);
a=loglog(x,px1,'color',[0 0 0.80392]);
hold on;
plot(x,px2);
b=loglog(x,px2,'color',[1 0.18824 0.18824]);
semilogx(x,px1);
set(gca,'XDir','reverse');
% set(gca,'linewidth',1,'xlim',[1,48],'xtick',1:1:48,'yminorgrid','off','tickdir','out','box','on','layer','top', 'FontName','Times New Roman','fontsize',16);
xlabel('Period/h','FontName','Times New Roman','fontsize',16);
ylabel('Power Spectral Density','FontName','Times New Roman','fontsize',16);
legend([a,b],'Zonal','Meridional','location','best');
% text(5.6,0.4,'d','color','r','fontname','Times New Roman','fontsize',16);
grid on;
%%
% Find the inflection point
j=gradient(px1,x); %Derivation
a=find(abs(j(1,1:floor(1025/2)))>0.1);
a1=find(abs(j(1,floor(1025/2):end))>0.1)+floor(1025/2)-1;
b=px1(a(end));
b1=px1(a1(1));
c1=x(a1(1));
%%
Pw=abs(fft(xn)).^2/500;
psd=abs(fft(xn)).^2/500;
%%
p=psd(xn);
评论0