用fft做。R里面有现成的FFT的函数,不需要另外写。
下面是用太阳黑子的数据做的频谱分析,原先是MATLAB的代码,我改成了R的代码,结果是R的频谱分析结果模型要比MATLAB的差,而且还有异常点要手工处理。
<br />
data(sunspot.year)<br />
s=sunspot.year</p>
<p>fs=1 #采样频率<br />
nyquist=fs/2</p>
<p>y=fft(s)<br />
n=length(y)<br />
power=abs(y[1:ceiling(n/2)])^2</p>
<p>freq=(1:ceiling(n/2))/ceiling(n/2)*nyquist<br />
plot(freq,power)</p>
<p>period=1/freq<br />
plot(period,power)<br />
</p>