フーリエ級数展開を可視化するMATLABスクリプトを書いてみた

何となくフーリエ級数展開を説明する可視化が欲しいなー，と思い立ちMATLABで書いてみました．類似の実装は世の中に山のようにあるはずですがあえて調べません．

出力はこんな感じになります．ライブスクリプトを実行するとアニメーションになってくれるので，MATLABをお持ちの方は下記ソースコードをコピペして実行してみてください．

clear
clc
close all

n=1:7; %項数の設定
b=zeros(max(n),1); %係数ベクトル
t=linspace(0,2*pi,200); %1周期の時刻

% fの設定
plot_f=[sign(sin(t))]; %矩形波
% plot_f=[t/pi];plot_f=plot_f-mean(plot_f); %のこぎり波
% plot_f=cumsum([sign(cos(t))])*(t(2)-t(1));
% plot_f=sin(t)+1/3*sin(3*t)+0.1*randn(size(t));
plot_f=plot_f-mean(plot_f); %平均を0にシフト

% 係数の数値積分(関数の内積)
for n1=n
    b(n1)=sin(n1*t)*plot_f'/(sin(n1*t)*sin(n1*t)');
end

figure
plot_t=[];plot_y=[];plot_x=[];
for t1=t
    xc=0;yc=0;
    clf
    axis equal;grid on;
    set(gca,'fontname','Arial','fontsize',12)

    hold on;
    for n1=n
        plot([xc,b(n1)*cos(n1*t1)+xc],[yc,b(n1)*sin(n1*t1)+yc],'b-o')
        plot(b(n1)*cos(n1*t)+xc,b(n1)*sin(n1*t)+yc,'k:')
        xc=xc+b(n1)*cos(n1*t1);
        yc=yc+b(n1)*sin(n1*t1);
    end
    plot_t=[plot_t;t1];
    plot_x=[plot_x;xc];
    plot_y=[plot_y;yc];
    plot(t, plot_f,'r--')
    plot([xc,plot_t(end)],[yc,yc],'b:o')
    plot(plot_x, plot_y,'b--')
    plot(plot_t, plot_y,'g-')
    hold off
    drawnow
end

設定するべき箇所は項数nと被展開関数の数値列plot_fです．時刻のベクトルtを1周期として$\sin nt$で展開したものを図示します．

展開したい関数としてplot_f=sin(t)+1/3*sin(3*t)+0.1*randn(size(t))なんかを選ぶと基本波と3倍高調波のみ現れて雑音が消去されていることなんかもわかります．

動画ファイルとしての書き出しはVideoWriterを参照してください．

aviをgifに変換したい場合，下記のプログラムが便利です．対話型で細かな設定(ループ回数，最初と最後の静止時間)が可能です．