3

MATLAB を使用して波形のフーリエ係数を計算しようとしています。係数は、次の式を使用して計算できます。

ここに画像の説明を入力

ここに画像の説明を入力

T は、オメガ = 2pi を与える 1 になるように選択されます。

ただし、積分の実行に問題があります。関数は、三角波 (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

あなたの問題は、 Mathworksのドキュメントsawtoothに次のように返されるものであるように見えます：

sawtooth(t,width) は、幅 (0 から 1 の間のスカラー パラメーター) が最大値が発生する 0 から 2π の間の点を決定する、修正された三角波を生成します。この関数は、間隔 0 から 2π の幅で -1 から 1 に増加し、間隔 2π から 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お役に立てれば！

于 2013-02-22T17:07:15.363 に答える
0

コードのトラブルシューティングを行うために、使用している関数をプロットして調査し、クワッド関数がそれらをどのようにサンプリングするかを調べます。それらをアンダーサンプリングしている可能性があるため、最小ステップ サイズが関数の周期より少なくとも 10 倍小さいことを確認してください。

Matlab に組み込まれている FFT を使用することをお勧めします。FFT はスペクトルを計算する最も効率的な方法であるだけでなく (n*log(n) は配列の長さ n に依存し、積分は n^2 に依存します)、周波数ポイントも自動的に得られます。 （等間隔の）時間データによってサポートされています。積分を自分で計算する場合 (データポイントが等間隔でない場合に必要になる場合があります)、解決されない周波数データを計算する可能性があります (間隔が 1 より狭い/時間の間隔を超える、つまり「フーリエ限界」を超えている)。

于 2013-02-22T16:37:43.160 に答える