3

我正在尝试使用 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 系数是针对三角波给出的，当 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
4

2 回答 2

3

在我看来，您的问题是sawtooth返回mathworks 文档的原因：

sawtooth(t,width) 生成修改后的三角波，其中宽度是 0 到 1 之间的标量参数，决定了 0 到 2π 之间出现最大值的点。该函数在 0 到 2π 宽度的区间上从 -1 增加到 1 ，然后在 2π 宽度到 2π 的区间上从 1 到 -1 线性减小。因此，参数 0.5 指定标准三角波，关于时间瞬间 π 对称，峰峰值幅度为 1。锯齿波 (t,1) 等效于锯齿波 (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. 我发现通过使用trapzwhich 效果很好，直到我使用的时间向量没有足够的分辨率，这让我相信这是一个数字问题而不是基本问题。希望这可以帮助！

于 2013-02-22T17:07:15.363 回答
0

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

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

于 2013-02-22T16:37:43.160 回答