PSD.m.zip_matlab 图像序列_pwelch_功率谱

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);