FFTの結果を複素数のまま立体グラフにしてみた

13399 ワード
フーリエ変換 FFT matlab matlab テキストリンク

第1章. やりたいこと

エクセルのデータ分析ツール、matlabのfft関数など、高速フーリエ変換の結果は複素数で得られることが多い。
これを絶対値に直し、共役複素を折り返すことで、周波数対パワー特性を得ることができる。
当然この方法では位相の情報は失われる。
それなら複素平面をの2次元と周波数の1次元を組み合わせた立体グラフをかけば、少々読みにくいが位相情報が失われない。(3Dプリンタで出したら面白いかも。)

第2章. やり方

周波数をxとする。複素フーリエ係数をy+ziとする。これでxyz空間を図示すればよい。

第3章. やってみた

例として、「横軸を数直線、縦軸を約数の個数」とするような離散値をFFTにかけてみよう。

3-1. 失敗例

次のコードを実行した。(218まで求めると3時間ほどかかる。210までだとすぐ終わる。またFFTでは、データの大きさは2の冪数である必要がある。)

ワークスペースに入力
n_max=2^18;
n=n_min:n_max;
result=divisorsCnts_progress(n,5);

sjze=size(result);
x_min=1;
x_max=sjze(1)*sjze(2);
x=x_min:x_max;
fs=1;%サンプリング周波数
T=1/fs;%サンプリング周期
L=x_max-x_min+1;
Y=fft(result);
P2=abs(Y/L);
P1=P2(1:L/2+1);
P1(2:end-1)=2*P1(2:end-1);
f=fs*(0:(L/2))/L;
plot(f,P1);
title('約数の個数の周期性(FFT)');
xlabel('f[/自然数1増加]');
ylabel('|P1(f)|');

f=fs*(0:L)/L;
f=f(1:end-1);
Y_=real(Y);
Z_=imag(Y);
plot3(f, Y_, Z_);
title('約数の個数の周期性及び位相(FFT)');
xlabel('f[/自然数1増加]');
ylabel('Re(P1)');
zlabel('Im(P1)');
divisorsCnt.m
function cnt = divisorsCnt( natural )
%DIVISORSCNT 入力された自然数の約数の個数を返す。
    cnt = size(divisors(natural));
    cnt = cnt(1)*cnt(2);
end
divisorsCnts_progress.m
function cnts = divisorsCnts_progress(naturals, progViewFlu)
    sjze = size(naturals);
    sjze = sjze(1)*sjze(2);
    fprintf("自然数の個数:%d",sjze);
    cnts(sjze)=0;
    tic;
    j=1;

    for i=1:sjze
        cnts(i)=divisorsCnt(i);
        if(toc>progViewFlu*j)
            fprintf("自然数%dの時点での時間:%d秒", naturals(i), toc);
            j=j+1;
        end

    end
end

これを観測すると、次のいくつかの図のように実部-f特性は半分の周波数に線対称、虚部-f特性は半分の周波数で値0の点に点対称となっている。

実は、半分以降は「負の周波数成分」のエイリアスなので、無効だ。
したがって、次のようにする。

3-2. 成功例

ワークスペースに入力
f=fs*(0:(L/2))/L;
Y_=real(Y);
Z_=imag(Y);
Y_=Y_(1:uint64(end/2)+1);
Z_=Z_(1:uint64(end/2)+1);
plot3(f, Y_, Z_);
title('約数の個数の周期性及び位相(FFT)');
xlabel('f[/自然数1増加]');
ylabel('Re(P1)');
zlabel('Im(P1)');

結果、次のようなグラフが得られた。

参考