Обратное преобразование Фурье

Question

Я пытаюсь написать алгоритм 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 имеет функцию Иффт(), но я хотел бы написать мой собственный)? Я использую для него этот код:

Результат должен выглядит как (EN ошибка):

Где черный участок и синий х^Участок х. Мой результат - прямая.

Answer 1

Я не знаю, где вы получили эту формулу в изображении с, но для N-точечного преобразования Фурье дискретного сигнала, вам нужно только просуммировать k от 0 к N-1, чтобы получить точную реконструкцию. N Функции ортогонального базиса могут восстанавливать любой N -мерный сигнал. Вы можете ввести в заблуждение DTFT с DFT (вы хотите, чтобы второй из них).

Я также не понимаю формулы, которые вы использовали для коэффициентов 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)