0

私はIFTアルゴリズムを書こうとしています。これが私のコードです：

%% Fourier Analysis
N = 20;
T1 = 0.25*N;
T0 = N;
x = [zeros(T1,1); ones(T0/2,1); zeros(T1,1)]';

omega = 2*pi*1/T0;

ak = zeros(1,2*N+1);

for k = -N:N
    if k == 0
        ak(k+N+1) = 2*T1/T0;
    else
        ak(k+N+1) = sin(k*omega*T1)/(k*pi);
    end
end

%Approximate the periodic symmetric square wave
t=linspace(-0.5,0.5,N);

for n=1:length(t)
    xN(n)=0;
    for k = -N:N
        xN(n) = xN(n) + ak(k+N+1).*exp(1i*k*omega*t(n));
    end
end

何が問題なのですか (Matlab には ifft() 関数があることは知っていますが、自分で書きたいと思います)? 私はそれにこのコードを使用します：

場合

結果は次のようになります (EN はエラーです)。

ここに画像の説明を入力

黒いプロットは xN で、青いプロットは x です。私の結果は直線です。

4

1 に答える 1

3

画像のどこからその式を取得したかはわかりませんが、離散信号の点フーリエ変換の場合、正確な再構成を得るためにからまでNを合計するだけで済みます。直交基底関数は、任意の次元の信号を再構築できます。DTFTとDFTを混同している可能性があります(これらの 2 番目が必要です)。k0N-1NN

akまた、係数に使用した式もわかりません。それらを使用して計算しsin、正弦波ではなく複雑な指数関数で合計します。

複雑な指数関数で離散フーリエ変換 (DFT) を実行する場合、コードは次のようになります。時間信号と複素基底関数ckの内積から係数を取得します。x

clear; clc;

N = 20;                                        %// number of points

x = [zeros(1,N/4),ones(1,N/2),zeros(1,N/4)];   %//time signal data
n = 0:N-1;

ck = zeros(1,N);
for k = 0:N-1                                  %//cacluate coefficients
    ck(k+1) = sum(x.*exp(-1i*2*pi*k*n/N));
end

xN = zeros(1,N);
for k = 0:N-1                                  %// add partial frequencies
    xN = xN + ck(k+1)*exp(1i*2*pi*k*n/N)/N;
end

plot(n,xN)

ここに画像の説明を入力

実正弦波を使用して通常のフーリエ級数を実行する場合、コードは次のようになります。

clear; clc;

N = 20;                                         %// number of points
T = 1;                                          %// fundamental frequency
omega = 2*pi/T;                                 %// angular frequency

t = linspace(-0.5,0.5,N);                       %// time points
x = [zeros(1,N/4),ones(1,N/2),zeros(1,N/4)];    %// time signal data

a0 = sum(x)/N;  ak = zeros(1,N);    bk = zeros(1,N);
for k = 1:N-1                                   %// calculate coefficients
    ak(k) = sum(x.*cos(omega*k*t))/N;
    bk(k) = sum(x.*sin(omega*k*t))/N;
end

xN = a0*cos(omega*0*t);                         %// add DC offset
for k = 1:N-1                                   %// add partial frequencies
    xN = xN + ak(k)*cos(omega*k*t) + bk(k)*sin(omega*k*t);
end

plot(t,xN)
于 2015-03-19T18:59:46.473 に答える