periodogram加显著性检验

伽蓝鸟

MATLAB的函数periodogram可以算出傅里叶功率谱，但是这个函数本身没法画显著性检验的线。看了下Wilks那本统计书，只要加一行命令就能标出显著性了。代码如下：

clc,clear,close all;fclose all;
%%  读入数据
fid=fopen('haikou2010.dat');
data=fread(fid,'single');fclose(fid);
%%  计算傅里叶谱及显著性
data=zscore(data);                                  % % 我都是用的标准化变量，请标准化
[Pxx,f]=periodogram(data,[],[],1,'power');          % % 频率请根据数据做修改
alpha=autocorr(data);alpha=alpha(2);                % % 滞后1时刻自相关，alpha=0则为白噪音
rednoise95=5.991*(1-alpha^2)./(length(data)*(1+alpha^2-2*alpha*cos(2*pi*f)));
                        % % 显著性检验就这一行代码
                        % % 注意如果上面的频率不是1，最后的cos(2*pi*f)要做相应修改
                        % % 卡方右尾分位数 95%对应着5.991，可以换成4.605(90%)或者9.210(99%)
%%  画图
figure,hold on;box on;
plot(1./f,Pxx,'-k','linewidth',2);
plot(1./f,rednoise95,'--r','linewidth',1.5);
set(gca,'xscale','log','xtick',[1 2 10 20 100],'xdir','reverse',...
    'fontsize',16,'ytick',0:0.1:0.4);
xlabel('Period (day)','fontweight','bold');
ylabel('Power','fontweight','bold');
title('Power Spectrum','fontweight','bold');
axis([1 100 0 0.4])

乔老师教的我们气候统计，就找了她的一篇文章作为验证（Qiao Y, Zhang C, Jian M. Role of the 10–20-day oscillation in sustained rainstorms over Hainan, China in October 2010[J]. Advances in Atmospheric Sciences, 2015, 32(3): 363-374.）。对比结果如下


另外，Aires前辈早就写过一个帖子讲功率谱(http://bbs.06climate.com/forum.php?mod=viewthread&tid=14978)。我也对比过了，结果很接近。

bbs

学习了！！数据标准化，怎么样才算标准化呢？

伽蓝鸟

bbs 发表于 2020-2-21 16:30
学习了！！数据标准化，怎么样才算标准化呢？

zscore就是标准化了啊。减去均值，除以标准差。

bbs

多谢多谢！！

xxr

您好，请问这个fs=1怎么理解呢，查资料看说是采样频率，对于咱们气象方面的，如果是日数据，这个fs是=1吗