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

关于傅里叶变换有几个cmets，希望能给大家解释一下。另外，我不知道您所说的“数学傅里叶变换”是什么意思，因为您的代码中没有一个表达式类似于 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 的每一项的复系数。请注意，它使用 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]')

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

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

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

警告：复数 X 和/或 Y 参数的虚部被忽略

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

norm(imag(x_approx))

 ans =

      1.402648396024229e-12

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