逆傅里叶变换
Inverse fourier transform
我正在尝试编写 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。我的结果是直线。
我不知道你从哪里得到图像中的公式,但是对于离散信号的 N
点傅立叶变换,你只需要从 k
求和 0
到 N-1
以获得精确的重建。 N
正交基函数可以重构任何 N
维信号。您可能会混淆 DTFT with the DFT(您想要其中的第二个)。
我也不明白您用于 ak
系数的公式。您用 sin
计算它们,然后用复指数而不是正弦波对它们求和。
如果你想用复指数做离散傅里叶变换 (DFT),代码应该类似于我下面的代码。您可以从时间信号 x
和复基函数的内积中得到 ck
系数。
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)
我正在尝试编写 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。我的结果是直线。
我不知道你从哪里得到图像中的公式,但是对于离散信号的 N
点傅立叶变换,你只需要从 k
求和 0
到 N-1
以获得精确的重建。 N
正交基函数可以重构任何 N
维信号。您可能会混淆 DTFT with the DFT(您想要其中的第二个)。
我也不明白您用于 ak
系数的公式。您用 sin
计算它们,然后用复指数而不是正弦波对它们求和。
如果你想用复指数做离散傅里叶变换 (DFT),代码应该类似于我下面的代码。您可以从时间信号 x
和复基函数的内积中得到 ck
系数。
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)