MATLAB|[MATLAB学习笔记]采用快速傅里叶变换求时间序列的周期项时间序列|快速傅里叶变换|ma

[MATLAB学习笔记]采用快速傅里叶变换求时间序列的周期项
1. 背景现有长度为11年的5个时间序列，为某拟研究对象的5个参数。现计划先通过快速傅里叶变换求每个系数序列的显著周期项，再分别按照傅里叶序模型列进行拟合，得到系数和时间的模型，即达到一个效果：输入一个时间，可以得到对应的5个参数。
2. 实现 2.1 序列中含有缺失天数，即为非等间隔采样 【MATLAB|[MATLAB学习笔记]采用快速傅里叶变换求时间序列的周期项】需先对序列进行插值计算，不同插值方法对结果的影响这里先不讨论。先用matlab中自带的spline函数对数据进行三次样条插值，得到完整序列。
（实际上，本文中采用的数据点总长度为4018天，缺失的天数为25天，缺失率为0.6%，认为不同插值方法的影响可以忽略）
插值部分代码如下：

clear%读取参数序列 load ./coefFile/nml_PolyCurveCoef_1998.182-2009.182.txt data = https://www.it610.com/article/nml_PolyCurveCoef_1998_182_2009_182; %总天数：s:50995(1998.7.1)e:55196(2009.12.31) 一共4201天%截取1998.7.1-2009.7.1之间的数据 %总天数：s:50995(1998.7.1)e:55013(2009.7.1) 一共4018天 data_temp = ones(4018,9); for i = 1: length(data) if data(i,1) == 2009 && data(i,2) == 7 && data(i,3) == 1 data_temp(i,:) = data(i,:); break; else data_temp(i,:) = data(i,:); end enddata = data_temp; it = 1; for i = 1: length(data) if data(i,5) == 0 && data(i,6) == 0 && data(i,7) == 0 is_zero_ind(it) = i; %is_zero_ind数组储存拟合参数为0的天数在矩阵中的行位置（即缺失数据天的位置） it = it + 1; endend mjd_time_temp = data(:,4); mjd_time_interp = mjd_time_temp(is_zero_ind); %读取需要插值的mjd mjd_time_interp = mjd_time_interp - 50965; %减去初始天数time_temp = data(:,1:3); time_interp = time_temp(is_zero_ind,1:3); %读取需要插值的年月日时间信息data1=data; [row,col]=size(data1); for j=1:length(is_zero_ind) data1(is_zero_ind(j),:)=ones(1,col)*NaN; end data1(isnan(data1))=[]; %此处删除完所有NaN后变成一列 data1=reshape(data1,length(data1)/col,col); %再按照固定的列数将其重新排列就可以得到新矩阵%插值情况： %总天数：s:50995(1998.7.1)e:55013(2009.7.1) 一共4018天 %现有数据：3993天 %插值长度：25天 %data1不含缺失数据mjd_time = data1(:,4); mjd_mjd0 = mjd_time - 50965; a0 = data1(:,5); a1 = data1(:,6); a2 = data1(:,7); a3 = data1(:,8); a4 = data1(:,9); a0_interp = interp1(mjd_mjd0,a0,mjd_time_interp,'spline'); a1_interp = interp1(mjd_mjd0,a1,mjd_time_interp,'spline'); a2_interp = interp1(mjd_mjd0,a2,mjd_time_interp,'spline'); a3_interp = interp1(mjd_mjd0,a3,mjd_time_interp,'spline'); a4_interp = interp1(mjd_mjd0,a4,mjd_time_interp,'spline'); %插值得到的缺失天的数据 data_missingData = https://www.it610.com/article/[time_interp mjd_time_interp + 50965 a0_interp a1_interp a2_interp a3_interp a4_interp]; data_interp = data; for inp = 1: length(is_zero_ind) data_interp(is_zero_ind(inp),1:9) =data_missingData(inp,1:9); endnml_poly_coef_1998_182_2009_182_interp = data_interp; save ./matFile/nml_poly_coef_1998_182_2009_182_interp nml_poly_coef_1998_182_2009_182_interp % %绘制图 % % % subplot(5,1,1) % plot(mjd_mjd0,a0,'b.') % hold on % plot(mjd_time_interp,a0_interp,'r*') % ylabel('a0') % xlim([0 4020]) % % subplot(5,1,2) % plot(mjd_mjd0,a1,'b.') % hold on % plot(mjd_time_interp,a1_interp,'r*') % ylabel('a1') % xlim([0 4020]) % % subplot(5,1,3) % plot(mjd_mjd0,a2,'b.') % hold on % plot(mjd_time_interp,a2_interp,'r*') % ylabel('a2') % xlim([0 4020]) % % subplot(5,1,4) % plot(mjd_mjd0,a3,'b.') % hold on % plot(mjd_time_interp,a3_interp,'r*') % ylabel('a3') % xlim([0 4020]) % % subplot(5,1,5) % plot(mjd_mjd0,a4,'b.') % hold on % plot(mjd_time_interp,a4_interp,'r*') % xlabel('td = mjd - 50965') % xlim([0 4020]) % ylabel('a4')

输入文件中，数据缺失的天数5个参数用均用0填上了，所以is_zero_ind数组是为了找到这些数据缺失位置
2.2 快速傅里叶变换数据插值后得到完整的序列，进行快速傅里叶变换。由于没有学习过信号相关知识，所以对该方面了解并不深刻，建议和我一样的先去读一下参考博客2和3中的详细描述。主要在FFT计算后想以下面的形式展示：

文章图片

其将参数时间序列进行傅里叶变换后，画出周期和幅值。横坐标为days，而非一般信号分析里面频率。这里先给出我计算的程序，再记录一下我的理解：

clear load ./matFile/nml_poly_coef_1998_182_2009_182_interp.mat data = https://www.it610.com/article/nml_poly_coef_1998_182_2009_182_interp; year = data(:,1); mon = data(:,2); day = data(:,3); mjd_time = data(:,4); mjd_mjd0 = mjd_time - 50965; a0 = data(:,5); a1 = data(:,6); a2 = data(:,7); a3 = data(:,8); a4 = data(:,9); x = a1'; figure ; subplot(2,1,1); plot(mjd_mjd0,x,'.') xlim([0 4020]) xlabel('td = mjd - 50965') ylabel('a0')N = length(x); %采样点数 4018 个点 time_len = mjd_time(end) - mjd_time(1); %时间长度（单位：天）4018天 T = time_len/N; % 采样周期等于采样长度除去采样点数(单位：天)表示每隔T天采集一个数据 Fs=1/T; % 采样频率（单位：每天）n = 0:N-1; Y = fft(x); % 计算出来的fft是个复数df=Fs/N; %频率间隔 f=(0:1:N/2)*df; % 频率刻度 period = 1./f; YY=Y(1:N/2+1); % 提取正频率的部分 YY(2:end-1)=2*YY(2:end-1); % 负频率合并到正频A=abs(YY); % 计算频域序列YY的幅值 Pha=angle(YY); % 计算频域序列YY的相角弧度subplot(2,1,2) plot(period,A)% 绘制频域序列YY的幅频图 grid on xlim([0 4020]) xlabel('period(days)') ylabel('Amplitude')

在一般信号分析的教程中，频谱图的横坐标为频率，单位为Hz（次/秒），本文中采用的时间序列，时间长度为4018天，采样点数为4018个，则采样间隔为1，即每1天采集一个数据，那么此时根据博客教程所示画出来的频谱图横坐标的“频率”单位就不是Hz（次/秒）了，而是（次/天），这个单位对应的”频率“在本例中是个非常小的数字，此时为了更好的展示时间序列中的周期项，只需要将横坐标的数字求倒数就可以了，频率的倒数是周期。则画出来的a1参数的周期-幅值图为：