MATLAB fft 与数学傅里叶

Question

我正在使用以下代码生成信号的 fft 和数学傅里叶变换。然后我想在数学上重新创建 fft 的原始信号。这适用于数学信号，但不适用于 fft，因为它是离散变换。有谁知道我可以对我的逆变换方程进行哪些更改，使其适用于 fft？

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');```

Answer 1

关于傅里叶变换的评论有几点，希望我能为你解释一切。另外，我不知道你所说的 "Mathematical Fourier transform" 是什么意思，因为你代码中表达式的 none 类似于 Fourier series of the sawtooth wave.

要准确理解fft function的作用，我们可以一步一步来。首先，按照您的代码，我们创建并绘制锯齿波的一个周期。

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');

我们现在可以计算一些东西了。首先，Nyquist frequency，样本的最大可检测频率。

f_max = 0.5*fs

 f_max =

    512

另外，最小可检测频率，

f_min = 1/t(end)

 f_min =

    1.000977517106549

现在用 MATLAB 函数计算离散傅里叶变换：

X = fft(x)/n;

此函数获取discrete Fourier transform的每一项的复数系数。请注意，它使用指数表示法计算系数，而不是根据正弦和余弦。除以n是为了保证第一个系数等于样本的算术平均值

如果要绘制转换信号的 magnitude/phase，可以键入：

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]')

这个情节是什么意思？它在图中显示了表示原始信号（锯齿波）的傅立叶级数中各项的复系数的幅度和相位。您可以使用此系数根据（截断的）傅立叶级数获得信号近似值。当然，要做到这一点，我们需要整个变换（不仅是前半部分，因为通常会绘制它）。

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 循环更好的方式完成，并且代码中会弹出一些警告，警告用户为速度预分配每个变量。

上面的代码使用截断傅里叶级数的项计算每次 ti 的 x(ti)。如果我们绘制原始信号和近似信号，我们得到：

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

原始信号和近似信号几乎相等。事实上，

norm(x-x_approx)

 ans =

      1.997566360514140e-12

几乎为零，但不完全为零。此外，由于在计算近似信号时使用了复系数，上图会发出警告：

Warning: Imaginary parts of complex X and/or Y arguments ignored

但是你可以检查虚数项非常接近于零。由于计算中的舍入误差，它不完全为零。

norm(imag(x_approx))

 ans =

      1.402648396024229e-12

请注意上面代码中如何解释和使用 fft 函数的结果以及它们如何以 exp 形式表示，而不是像您编码的那样用正弦和余弦表示。

MATLAB fft 与数学傅里叶

MATLAB fft vs Mathematical fourier

matlab

fft