基于FFT的Matlab离散时间序列正弦余弦和拟合与预测技术问询

核心思路 #

离散实值时间序列可分解为不同频率的正弦、余弦项加权和，FFT能直接给出这些频率对应的系数——基于最小平方误差准则，FFT的结果就是最优拟合的系数。无需额外解线性方程组（若手动构建矩阵求解，结果与FFT完全一致）。

代码问题修正 #

你当前的代码存在几个关键问题：

• 变量名冲突：用x同时存储时间点和求解结果，导致数据覆盖
• 设计矩阵A构建错误：对角矩阵无法正确表示各频率项对每个时间点的贡献，正确的A应该是每行对应一个时间点，每列对应一个余弦/正弦基函数
• 频率与时间的对应关系错误：未正确使用采样时间间隔计算相位

正确实现方案 #

方案1：直接用FFT重构拟合序列（最简便） #

FFT本身就是对离散傅里叶级数的最优拟合，直接利用FFT结果重构即可：

% 加载数据（假设eta_train是时间序列，t_train是对应的时间向量）
b = eta_train;
t = t_train;
N = length(b);

% 计算采样频率和时间间隔
Delta_t = t(2) - t(1); % 时间间隔
Fs = 1 / Delta_t;      % 采样频率

% 计算FFT及对应频率
B = fft(b);
frequencies = (0:N-1) * Fs / N;

% 计算余弦项系数alpha和正弦项系数beta
alpha = 2 * real(B) / N;
beta = -2 * imag(B) / N;
% 直流分量（k=0）和Nyquist频率（k=N/2，当N为偶数时）不需要乘2
alpha(1) = real(B(1))/N;
if mod(N,2) == 0
    alpha(N/2 + 1) = real(B(N/2 + 1))/N;
end

% 重构拟合序列
b_fit = zeros(size(b));
for k = 1:N
    b_fit = b_fit + alpha(k)*cos(2*pi*frequencies(k)*t) + beta(k)*sin(2*pi*frequencies(k)*t);
end

% 绘制原始序列和拟合序列
figure;
plot(t, b, 'b-', t, b_fit, 'r--');
legend('原始序列', '拟合序列');
xlabel('时间');
ylabel('序列值');

方案2：手动构建设计矩阵求解线性方程组 #

若想手动构建矩阵A并求解系数，本质和FFT结果一致，代码如下：

b = eta_train;
t = t_train;
N = length(b);
Delta_t = t(2)-t(1);
Fs = 1/Delta_t;
frequencies = (0:N-1)*Fs/N;

% 构建设计矩阵A：每行是一个时间点的余弦、正弦项值
% 注意：对于实序列，频率是共轭对称的，只需取前半部分即可减少计算量
num_freq = ceil(N/2);
A = zeros(N, 2*num_freq - (mod(N,2)==0)); % 直流分量只有余弦项

% 直流分量（k=0）
A(:,1) = cos(2*pi*frequencies(1)*t);

% 其他频率分量
col_idx = 2;
for k = 2:num_freq
    A(:,col_idx) = cos(2*pi*frequencies(k)*t);
    A(:,col_idx+1) = sin(2*pi*frequencies(k)*t);
    col_idx = col_idx + 2;
end

% 求解系数x
x = A \ b;

% 重构拟合序列
b_fit = A * x;

% 绘图对比
figure;
plot(t, b, 'b-', t, b_fit, 'r--');
legend('原始序列', '拟合序列');
xlabel('时间');
ylabel('序列值');

关键说明 #

• 对于实值序列，FFT结果是共轭对称的，因此只需保留前半部分频率即可，后半部分冗余，能大幅减少计算量
• 两种方案得到的拟合结果完全一致，FFT方案效率更高，适合大序列
• 若要进行预测，只需将时间t扩展到未来时刻，代入拟合公式即可得到预测值

内容的提问来源于stack exchange，提问作者oviearies