如何使用 MATLAB 计算傅立叶系数

Posted 2023-03-31

【中文标题】如何使用 MATLAB 计算傅立叶系数

【英文标题】：How to compute Fourier coefficients with MATLAB

【发布时间】：2013-02-22 16:28:28

【问题描述】：

我正在尝试使用 MATLAB 计算波形的傅立叶系数。可以使用以下公式计算系数：

T 被选择为 1，这使得 omega = 2pi。

但是，我在执行积分时遇到问题。函数是三角波（如果我没记错的话，可以使用sawtooth(t,0.5) 生成）以及方波。

我已尝试使用以下代码（对于三角波）：

    function [ a0,am,bm ] = test( numTerms )
            b_m = zeros(1,numTerms);
            w=2*pi;
            for i = 1:numTerms
                    f1 = @(t) sawtooth(t,0.5).*cos(i*w*t);
                    f2 = @(t) sawtooth(t,0.5).*sin(i*w*t);
                    am(i) = 2*quad(f1,0,1);
                    bm(i) = 2*quad(f2,0,1);
            end
    end

但是，它并没有接近我需要的值。 b_m 系数针对 a 三角波，当 m 为奇数时，从正项开始交替变化，应该是 1/m^2 和 -1/m^2。

对我来说主要问题是我不太了解积分在 MATLAB 中的工作原理，而且我不确定我选择的方法是否有效。

编辑： 澄清一下，这是我希望在确定系数时编写函数的形式：

这是使用 fft 的尝试：

    function [ a0,am,bm ] = test( numTerms )
            T=2*pi;
            w=1;
            t = [0:0.1:2];
            f = fft(sawtooth(t,0.5));
            am = real(f);
            bm = imag(f);
            func = num2str(f(1));
            for i = 1:numTerms
                    func = strcat(func,'+',num2str(am(i)),'*cos(',num2str(i*w),'*t)','+',num2str(bm(i)),'*sin(',num2str(i*w),'*t)');
            end
            y = inline(func);
            plot(t,y(t));
    end

【问题讨论】：

fft() 有什么问题

我需要一个适合使用正弦和余弦近似函数的表格上的系数，我真的不知道如何从 fft 给出的答案中得到该表格，所以我认为这可能是一个更简单的方法。

FFT 在被积函数中使用表达式 exp(ix) = cos(x) + i sin(x)，因此要获得 cos 和 sin 部分，您只需要取实部和虚部即可.

编辑了一个包含 FFT 尝试的新示例，但仍无法按预期工作。

【参考方案1】：

在我看来，您的问题是 sawtooth 返回的 mathworks documentation 表示：

sawtooth(t,width) 生成修改后的三角波，其中宽度（0 到 1 之间的标量参数）确定 0 到 2π 之间出现最大值的点。该函数在 0 到 2πwidth 的区间上从 -1 增加到 1，然后在 2πwidth 到 2π 的区间上从 1 到 -1 线性减小。因此，参数 0.5 表示标准三角波，关于时间点 π 对称，峰峰值幅度为 1。sawtooth(t,1) 等价于 sawtooth(t)。

所以我猜这是您问题的部分。

在你回复之后，我又调查了一些。在我看来这是quad 函数；不是很准确！我把问题改写成这样：

    function [ a0,am,bm ] = sotest( t, numTerms )
      bm = zeros(1,numTerms);
      am = zeros(1,numTerms);
      % 2L = 1
      L = 0.5;
      for ii = 1:numTerms
        am(ii) = (1/L)*quadl(@(x) aCos(x,ii,L),0,2*L);
        bm(ii) = (1/L)*quadl(@(x) aSin(x,ii,L),0,2*L);
      end
      ii = 0;
      a0 = (1/L)*trapz( t, t.*cos((ii*pi*t)/L) );  
      % now let's test it
      y = ones(size(t))*(a0/2);
      for ii=1:numTerms
        y = y + am(ii)*cos(ii*2*pi*t);
        y = y + bm(ii)*sin(ii*2*pi*t);
      end
      figure; plot( t, y);
    end

    function a = aCos(t,n,L)
      a = t.*cos((n*pi*t)/L);
    end

    function b = aSin(t,n,L)
      b = t.*sin((n*pi*t)/L);
    end

然后我这样称呼它：

[ a0,am,bm ] = sotest( t, 100 );

我得到了：

甜蜜！！！

我真正改变的只是从 quad 到 quadl. 我通过使用 trapz 发现了这一点，直到我使用的时间向量没有足够的分辨率，这让我相信这是一个数字问题而不是基本问题。希望这会有所帮助！

【参考方案2】：

为了对您的代码进行故障排除，我将绘制您正在使用的函数并调查 quad 函数如何对它们进行采样。您可能对它们进行了欠采样，因此请确保您的最小步长比函数的周期小至少 10 倍。

我建议使用 Matlab 内置的 FFT。 FFT 不仅是计算频谱的最有效方法（它的 n*log(n) 取决于数组的长度 n，而 n^2 中的积分取决于），它还会自动为您提供频点您的（等间隔的）时间数据支持的。如果您自己计算积分（如果数据点的间距不相等，则可能需要），您可能会计算未解析的频率数据（间距小于 1/时间间距，即超出“傅立叶限制”）。