如何从已计算的DFT结果求解信号的正余弦形式傅里叶系数

前置修正：现有代码的小问题 #

注意：你之前代码里的n=0:A存在长度不匹配问题，A是x的长度16，所以n会生成0到16共17个元素，和16个元素的x不对应，需要修改为n=0:A-1 #

核心原理 #

对于长度为L的实序列x[n]（此处L=16），我们可以通过和信号等长的L点DFT结果，直接将信号拆解为直流分量+多个不同频率的正弦、余弦分量之和：

对于偶数长度L，分解公式为：
$$x[n] = \frac{X[0]}{L} + \frac{2}{L}\sum_{k=1}^{L/2-1}\left( Re(X[k])cos(\frac{2\pi kn}{L}) - Im(X[k])sin(\frac{2\pi kn}{L}) \right) + \frac{X[L/2]}{L}cos(\pi n)$$
其中X[k]是x[n]的L点DFT结果

具体实现步骤 #

• 第一步：计算16点DFT，你之前用的N=50是插值后的DFT，要分解原始信号必须使用和信号长度一致的N=16点DFT，才能保证分解无误差
• 第二步：按上述公式提取直流、各频率cos分量系数、sin分量系数、奈奎斯特频率分量系数
• 第三步：组合系数得到完整的分解表达式，也可以直接用系数重构信号验证准确性

可直接运行的Matlab代码 #

% 原始信号
x=[8.57837e-08, 2.07482e-06, 4.43796e-06, 7.66462e-06, 1.10232e-05, 1.35811e-05, 1.27958e-05, 5.94217e-06, 2.49168e-08, -6.58389e-06, -1.30551e-05, -1.345e-05, -1.07471e-05, -7.38637e-06, -4.42876e-06, -1.88811e-06 ];
L = length(x); % L=16
n = 0:L-1;

% 计算16点DFT
X = fft(x, L);

% 拆解为余弦+正弦分量
dc = X(1)/L; % Matlab索引从1开始，对应k=0的直流分量
fprintf('信号分解表达式：\n');
fprintf('x[n] = %.8e', dc);

% 存储各频率分量系数
cos_coef = zeros(1, L/2+1);
sin_coef = zeros(1, L/2+1);
cos_coef(1) = dc;

% 遍历中间频率分量k=1到7
for k = 1:L/2-1
    real_part = real(X(k+1));
    imag_part = imag(X(k+1));
    cos_k = 2*real_part/L;
    sin_k = -2*imag_part/L;
    cos_coef(k+1) = cos_k;
    sin_coef(k+1) = sin_k;
    if sin_k >=0
        fprintf(' + %.8e*cos(2π*%d*n/16) + %.8e*sin(2π*%d*n/16)', cos_k, k, sin_k, k);
    else
        fprintf(' + %.8e*cos(2π*%d*n/16) - %.8e*sin(2π*%d*n/16)', cos_k, k, abs(sin_k), k);
    end
end

% 添加奈奎斯特分量k=8，对应频率π
nyquist_coef = real(X(L/2+1))/L;
cos_coef(L/2+1) = nyquist_coef;
fprintf(' + %.8e*cos(π*n)\n', nyquist_coef);

% 验证分解准确性：重构信号和原始信号的最大误差
x_recon = zeros(1, L) + dc;
for k=1:L/2-1
    x_recon = x_recon + cos_coef(k+1)*cos(2*pi*k*n/L) + sin_coef(k+1)*sin(2*pi*k*n/L);
end
x_recon = x_recon + nyquist_coef*cos(pi*n);
max_error = max(abs(x - x_recon))

结果说明 #

• 运行代码后会直接打印出完整的正余弦组合表达式，你可以直接使用输出的系数
• 最后的max_error结果会在1e-16量级，属于浮点计算误差，说明分解和原始信号完全一致
• 如果你需要振幅-相位形式的表达式，可以将每组cos、sin分量合并为A_k*cos(2πkn/16 + φ_k)，其中A_k = sqrt(cos_coef(k)^2 + sin_coef(k)^2)，φ_k = atan2(-sin_coef(k), cos_coef(k))

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