fft - 如何使用 MATLAB 计算傅里叶系数

Question

我正在尝试使用 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

Answer 1

在我看来，您的问题是 sawtooth 返回 mathworks documentation说:

sawtooth(t,width) generates a modified triangle wave where width, a scalar parameter between 0 and 1, determines the point between 0 and 2π at which the maximum occurs. The function increases from -1 to 1 on the interval 0 to 2πwidth, then decreases linearly from 1 to -1 on the interval 2πwidth to 2π. Thus a parameter of 0.5 specifies a standard triangle wave, symmetric about time instant π with peak-to-peak amplitude of 1. sawtooth(t,1) is equivalent to 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 找到了它使用没有足够的分辨率，这让我相信这是一个数字问题而不是基本问题。希望这对您有所帮助!