MATLAB FFT против математического Фурье - PullRequest
Новая
       10

MATLAB FFT против математического Фурье

0 голосов
/ 07 октября 2019

Я использую следующий код для генерации БПФ и математического преобразования Фурье сигнала. Я хочу затем математически воссоздать оригинальный сигнал БПФ. Это работает с математическим сигналом, но не с БПФ, поскольку это дискретное преобразование. Кто-нибудь знает, какое изменение я могу внести в мое уравнение обратного преобразования, которое заставит его работать на БПФ? 

clear all; clc;
N = 1024;
N2 = 1023;
SNR = -10;
fs = 1024;
Ts = 1/fs;
t = (0:(N-1))*Ts;
x = 0.5*sawtooth(2*2*pi*t);
x1 = fft(x); 
Magnitude1 = abs(x1);
Phase1 = angle(x1)*360/(2*pi);

for m = 1:1024
   f(m) = m;                % Sinusoidal frequencies
   a = (2/N)*sum(x.*cos(2*pi*f(m)*t));      % Cosine coeff. 
   b = (2/N)*sum(x.*sin(2*pi*f(m)*t));       % Sine coeff
   Magnitude(m) = sqrt(a^2 + b^2);                % Magnitude spectrum
   Phase(m) = -atan2(b,a);                   % Phase spectrum
end

subplot(2,1,1);
plot(f,Magnitude1./512);   % Plot magnitude spectrum
       ......Labels and title.......
subplot(2,1,2);
plot(f,Magnitude,'k');    % Plot phase spectrum
ylabel('Phase (deg)','FontSize',14);
pause();

x2 = zeros(1,1024);         % Waveform vector 
for m = 1:24
    f(m) = m;               % Sinusoidal frequencies
    x2 = (1/m)*(x2 + Magnitude1(m)*cos(2*pi*f(m)*t + Phase1(m)));  
end
x3 = zeros(1,1024);         % Waveform vector
for m = 1:24
    f(m) = m;               % Sinusoidal frequencies
    x3 = (x3 + Magnitude(m)*cos(2*pi*f(m)*t + Phase(m)));  
end
plot(t,x,'--k'); hold on;
plot(t,x2,'k');
plot(t,x3,'b');```

1 Ответ

0 голосов
/ 07 октября 2019

Есть несколько комментариев о преобразовании Фурье, и я надеюсь, что смогу объяснить все для вас. Кроме того, я не знаю, что вы подразумеваете под «математическим преобразованием Фурье», так как ни одно из выражений в вашем коде не похоже на ряд Фурье пилообразной волны .

Чтобы точно понятьчто делает fft функция , мы можем делать шаг за шагом. Сначала, следуя вашему коду, мы создаем и строим один период пилообразной волны. 

n = 1024;
fs = 1024;
dt = 1/fs;
t = (0:(n-1))*dt;
x = 0.5*sawtooth(2*pi*t);

figure; plot(t,x); xlabel('t [s]'); ylabel('x');

enter image description here

Теперь мы можем вычислить несколько вещей. Во-первых, частота Найквиста , максимальная обнаруживаемая частота из выборок. 

f_max = 0.5*fs

 f_max =

    512

Также минимальная обнаруживаемая частота 

f_min = 1/t(end)

 f_min =

    1.000977517106549

Рассчитайте теперь дискретное преобразование Фурье с функцией MATLAB: 

X = fft(x)/n;

Эта функция получает комплексные коэффициенты каждого члена дискретного преобразования Фурье . Обратите внимание, что он рассчитывает коэффициенты, используя выражение exp, а не в виде синусов и косинусов. Деление на n должно гарантировать, что первый коэффициент равен среднему арифметическому для выборок

. Если вы хотите построить амплитуду / фазу преобразованного сигнала, вы можете набрать: 

f = linspace(f_min,f_max,n/2); % frequency vector
a0 = X(1); % constant amplitude
X(1)=[]; % we don't have to plot the first component, as it is the constant amplitude term
XP = X(1:n/2); % we get only the first half of the array, as the second half is the reflection along the y-axis

figure
subplot(2,1,1)
plot(f,abs(XP)); ylabel('Amplitude');
subplot(2,1,2)
plot(f,angle(XP)); ylabel('Phase');
xlabel('Frequency [Hz]')

enter image description here

Что означает этот график? На рисунке показаны амплитуда и фаза комплексных коэффициентов членов в ряду Фурье, которые представляют исходный сигнал (пилообразная волна). Вы можете использовать эти коэффициенты для получения аппроксимации сигнала в терминах (усеченного) ряда Фурье. Конечно, для этого нам понадобится целое преобразование (а не только первая половина, как это обычно делается для построения графика). 

X = fft(x)/n;
amplitude = abs(X);
phase = angle(X);
f = fs*[(0:(n/2)-1)/n (-n/2:-1)/n]; % frequency vector with all components

% we calculate the value of x for each time step
for j=1:n
    x_approx(j) = 0;
    for k=1:n % summation done using a for
        x_approx(j) = x_approx(j)+X(k)*exp(2*pi*1i/n*(j-1)*(k-1));
    end
    x_approx(j) = x_approx(j);
end

Примечание: приведенный выше код для пояснения и делаетне намерены быть хорошо закодированы. Суммирование может быть выполнено в MATLAB гораздо лучше, чем использование цикла for, и в коде появятся некоторые предупреждения, предупреждающие пользователя о предварительном выделении каждой переменной для скорости.

Приведенный выше кодвычисляет x(ti) для каждого времени ti, используя условия усеченного ряда Фурье. Если мы построим и исходный сигнал, и аппроксимированный, мы получим: 

figure
plot(t,x,t,x_approx)
legend('original signal','signal from fft','location','best')

enter image description here

Исходный сигнал и приблизительный сигнал почти равны. На самом деле, 

norm(x-x_approx)

 ans =

      1.997566360514140e-12

Почти ноль, но не совсем ноль. Кроме того, на графике выше будет выдано предупреждение из-за использования комплексных коэффициентов при расчете приближенного сигнала:

Предупреждение: мнимые части комплексных аргументов X и / или Y игнорируются

Но вы можете проверить, что мнимый член очень близок к нулю. Это не совсем ноль из-за ошибок округления в вычислениях. 

norm(imag(x_approx))

 ans =

      1.402648396024229e-12

Обратите внимание на приведенные выше коды, как интерпретировать и использовать результаты функции fft и как онипредставлены в виде опыта, а не в виде синусов и косинусов, как вы кодировали.

...