我的功率谱图怎么这样？

阿吉木柯

图一

图二

我之前模仿学姐的基于傅立叶变换的功率谱图（图１），画了一个我的（图２）
我的为什么这样？
程序也是拷的，改了以下参数。：
附上程序：



clear all;
load jupin.mat;
x1 =AOzong;
x2 =Tzong;
x3 =Rzong;
n=732;
for j=1:12; by(j)=sum(AOzong(1:n/12,j))/n*12;end;%%多年月平均值
  for i=1:n/12; for j=1:12; Ryjp(i,j)=AOzong(i,j)-by(j); end; end;
R=Ryjp'; x1=R(:);
for j=1:12; by(j)=sum(Tzong(1:n/12,j))/n*12;end;%%多年月平均值
  for i=1:n/12; for j=1:12; Ryjp(i,j)=Tzong(i,j)-by(j); end; end;
R=Ryjp'; x2=R(:);
for j=1:12; by(j)=sum(Rzong(1:n/12,j))/n*12;end;%%多年月平均值
  for i=1:n/12; for j=1:12; Ryjp(i,j)=Rzong(i,j)-by(j); end; end;
  R=Ryjp'; x3=R(:);
%Fourier transform
%%距平序列资料标准化
x1 = x1./std(x1); x2 = x2./std(x2); x3 = x3./std(x3);

%% 功率谱估计及检验
%%Estimate the power spectral density (PSD) of a signal using a periodogram;
%%===[Pxx,f] = periodogram(x,window,nfft,fs);===
%The maxlag  m = n/2 ; k = 2*m.*f ; Tao = 2*m./k;
[Pxx,f] = periodogram(x1,hanning(n),n); %% same as the Pxx1; 可指定频率fs＝4
[Pxx1,k]=psd(x1,n,n/2,hanning(n),n/2,0.95,'none');
[Pxx2,k]=psd(x2,n,n/2,hanning(n),n/2,0.95,'none');
[Pxx3,k]=psd(x3,n,n/2,hanning(n),n/2,0.95,'none');
y1 = fft(x1,n); y2 = fft(x2,n); y3 = fft(x3,n);

norm1 = n/(2*std(x1)^2); y1(1) = []; norm2 = n/(2*std(x2)^2); y2(1) = [];
norm3 = n/(2*std(x3)^2); y3(1) = [];
power1 = abs(y1(1:n/2)).^2;  power2 = abs(y2(1:n/2)).^2;  %rows=nfft/2
power3 = abs(y3(1:n/2)).^2;  
Pyy1 = y1.* conj(y1) / n*2; Pyy2 = y2.* conj(y2) / n*2;
Pyy3 = y3.* conj(y3) / n*2;
nyquist = 1/2; freq = (1:n/2)/(n/2)*nyquist;
period = 1./f(2:n/2+1);  k1 = k(2:n/2+1); per = n/2./k1;
pow1 = power1./norm1; pow2 = power2./norm2;
pow3 = power3./norm3;

Alpha = 0.72;
for i = 1: n/2;
  Pk1(i) = sum(pow1)/n*(1-Alpha^2)/(1+Alpha^2-2*Alpha*cos(2*pi*k(i)/n));
  Pk2(i) = sum(pow2)/n*(1-Alpha^2)/(1+Alpha^2-2*Alpha*cos(2*pi*k(i)/n));
  Pk3(i) = sum(pow3)/n*(1-Alpha^2)/(1+Alpha^2-2*Alpha*cos(2*pi*k(i)/n));
  

end
Pk951 = Pk1.*((7.815+9.488)/2/3.5); Pk952 = Pk2.*((7.815+9.488)/2/3.5);% For m = N/2;
Pk953 = Pk3.*((7.815+9.488)/2/3.5);

%% 估计凝聚函数平方值(magnitude squared Coherence Function) between two signals
%%Estimate magnitude squared coherence function between two signals;
%%===Cxy = cohere(x,y,nfft,fs,window,numoverlap);===
[Cxy12,f] = cohere(x1,x2,n,2,hanning(n/2),0); %% fs 默认值为 2；
[Cxy13,f] = cohere(x1,x3,n,2,hanning(n/2),0);

%%window is a periodic Hanning window of length nfft,numoverlap plot(f,[Pxy Pxyc])= 0;
%%Pxyc is the testing of confidence Alpha=0.05.
%%===[Pxy,Pxyc,f] = csd(x,y,nfft,fs,window,numoverlap,p);===
[Gxy12,Pxyc12,f]=csd(x1,x2,n,2,hanning(n/2),0,0.95,'none'); %% 位相谱值随fs变化
[Gxy13,Pxyc13,f]=csd(x1,x3,n,2,hanning(n/2),0,0.95,'none');

%% 用时间长度L表示位相落后L=lag*sita/k/pi，单位为月。k = 2*m.*f
sita12 = atan(imag(Gxy12)./real(Gxy12));
sita13 = atan(imag(Gxy13)./real(Gxy13));
Ls12 = sita12(2:n/2+1)/2/pi./f(2:n/2+1);
Ls13 = sita13(2:n/2+1)/2/pi./f(2:n/2+1);

%% For Output
CSD = [real(Gxy12) real(Gxy13) ];
Pxyc = [real(Pxyc12) real(Pxyc13) ];
%%================================图形输出 =====================================
figure('color',[1 1 1])
subplot('position',[0.1 0.73 0.4 0.25])
Xticks = 2.^(fix(log2(min(period))):fix(log2(max(period))));
plot(period,pow1,'k-',period,Pk1,'k-.')%,period,Pk951,'-.')
xlabel('周 期 /a  （AO指数）')
set(gca,'YLim',[0,max(pow1)+1]);
%set(gca,'XLim',log2([min(Xticks),max(Xticks)]));
set(gca,'XLim',[log2(min(Xticks)+0.75),log2(max(Xticks))]);
%set(gca,'XDir','default');
set(gca,'XDir','reverse');
set(gca,'XTick',log2(Xticks(:)));
set(gca,'XTickLabel',Xticks);
ylabel('功 率 谱')
%title('Periodogram')
subplot('position',[0.1 0.40 0.4 0.25])
Xticks = 2.^(fix(log2(min(period))):fix(log2(max(period))));
plot(period,pow2,'k-',period,Pk2,'k-.')%,period,Pk952,'-.')
xlabel('周 期 /a  （平均气温）')
set(gca,'YLim',[0,max(pow2)-14]);
%set(gca,'XLim',log2([min(Xticks),max(Xticks)]));
set(gca,'XLim',[log2(min(Xticks)+0.75),log2(max(Xticks))]);
%set(gca,'XDir','default');
set(gca,'XDir','reverse');
set(gca,'XTick',log2(Xticks(:)));
set(gca,'XTickLabel',Xticks);
ylabel('功 率 谱')
%title('Periodogram')
subplot('position',[0.1 0.07 0.4 0.25])
Xticks = 2.^(fix(log2(min(period))):fix(log2(max(period))));
plot(period,pow3,'k-',period,Pk3,'k-.')%,period,Pk953,'-.')
xlabel('周 期 /a  （降水量）')
set(gca,'YLim',[0,max(pow3)]);
%set(gca,'XLim',log2([min(Xticks),max(Xticks)]));
set(gca,'XLim',[log2(min(Xticks)+0.75),log2(max(Xticks))]);
%set(gca,'XDir','default');
set(gca,'XDir','reverse');
set(gca,'XTick',log2(Xticks(:)));
set(gca,'XTickLabel',Xticks);
ylabel('功 率 谱')
%title('Periodogram')
break

阿吉木柯

最近被小波弄晕了呀！求指教！

njzqxt

楼主，请问原始资料是什么格式的呢？

lfys

支持分享精神

xueyanhou

楼主，你这个程序可以得到想要的结果了吗？

阿吉木柯

xueyanhou 发表于 2012-6-27 14:53

                               
登录/注册后可看大图

楼主，你这个程序可以得到想要的结果了吗？

没有，后来我就用的小波功率谱

yushaotian

阿吉木柯 发表于 2012-6-28 12:53

                               
登录/注册后可看大图

没有，后来我就用的小波功率谱

楼主，你好，我最近也在做功率谱，但是出来的图也奇奇怪怪的，所以麻烦你可不可以把你的小波功率谱也给我一份，谢谢哦~

斥鷃

LZ把横坐标拉开来试试~你学姐从32开始，你是从64开始的，可能有影响~还有你的标准化没有做好，应该要减去平均值的~不过感觉影响不是很大，看这段程序我服了，我那时候做的时候是按照课本手写程序上去的，结果有这些函数~~~

kongfeng0824

阿吉木柯 发表于 2012-6-28 12:53

                               
登录/注册后可看大图

没有，后来我就用的小波功率谱

可以贴一下小波功率谱的程序吗？

kongfeng0824

阿吉木柯 发表于 2012-5-15 14:24

                               
登录/注册后可看大图

最近被小波弄晕了呀！求指教！

楼主是分析降水和温度与AO、ENSO的关系吗？