如何用matlab正确绘制傅里叶级数的相位谱？

Question

我已经制作了一个基本周期为 2 秒的矩形脉冲 x 来进行傅立叶变换

t=-2:0.01:2; 
dt=t(2)-t(1); %increment of time
fs=1/dt % sampling rate 
n=length(t) %number of samples
X=fftshift(fft(x,n))/n %fourier transform of x%

f=linspace(-fs/2,fs/2, n) %making frequency axis

X_angle=angle(X); %the phase of X

我期待相位谱交替出现 -pi/2 和 pi/2，

但图表（太糟糕了，由于缺乏我的声誉，我无法发布它）

显示 X_angle 随着频率的增加而逐渐增加，范围从 -pi 到 pi

当我用

绘制幅度谱时，它工作得很好

plot(f,X_mag), X_mag=abs(X)

我想我必须对angle(X)进行一些幅度调整

但不是角度与 X 的大小无关吗？

我不知道为什么X_angle的绝对值随着绝对频率的增加而增加。

Answer 1

您在代码中犯了几个错误。

1. 您没有在代码中打印出 x 的值。 我假设你想分析一个对称的矩形脉冲。 但是，您的样本数量是奇数。 这意味着，您的信号不是对称的，这会导致小的 与理论教科书值相比存在差异。

   1. 矩形脉冲具有无限带宽。 但是，对于 FFT，采样率必须至少是最高频率的两倍 的信号。这意味着您的采样率必须至少为 2*infinity。 运气不好，你不能这样做。没有人能做到这一点。 结果，您将获得别名，这意味着您的结果包含错误。 好消息是，在矩形函数的情况下，这种影响可以得到补偿 在 sinc 函数的帮助下。

   2. 如果你做对了，你就会得到正确的教科书系数。 如果您的输入信号的形式类似于 (N=10) 11111-1-1-1-1-1，则函数是 奇函数。这意味着 f(-t) = -f(t)。在这种情况下，矩形可以由 一系列正弦函数。理论函数为：

   f(t) = 4/pi( sin(wt) + 1/3 sin(3wt) + 1/5 sin(5wt) + 1/7 sin(7wt) ...)

   没有像 sin(2wt) 或 sin(4wt) 这样的偶数频率。 这意味着，频谱中每隔一个频率的值为零。 由于数字噪声，这些值并不完全为零，而是接近于零。 从这些值计算相位会产生无意义的值。 其他频率是正弦函数的傅立叶变换值。 正弦函数的 FFT 是纯虚数。因为和的所有元素都有 符号相同，所有角度的值都相同，即 pi/2。

下面是修改后的代码：

close all;
clear all;
clc;

t=-2:0.01:(2-0.01); 
dt=t(2)-t(1); %increment of time
fs=1/dt; % sampling rate 
n=length(t); %number of samples
x = [ones(1,n/2), -ones(1,n/2)];
X=fft(x)/n; %fourier transform of x%

f=0:n-1; %making frequency axis
sincComp = @(N) (exp(-i*pi*(f/N)).*sinc(f/N)).';
%  Perform the compensation
comp = transpose(sincComp(n));
Xcorrected = X.*comp;

X_angle=angle(Xcorrected(2:2:n)); %the phase of X

figure;
subplot(3,1,1);
hold on;
stem(f(1:10), pi*abs(X(1:10))/2, 'ob', 'LineWidth', 3);
plot(f(1:10), pi*abs(Xcorrected(1:10))/2, 'or', 'LineWidth', 3);
hold off;
grid on;
title('Absolute value of FFT result', 'FontSize', 18);
xlabel('frequency', 'FontSize', 18);
ylabel('abs', 'FontSize', 18);
legend(['FFT'], ['FFT compensated']);
subplot(3,1,2);
hold on;
stem(f(1:10), pi*real(X(1:10))/2, 'ob', 'LineWidth', 3);
plot(f(1:10), pi*real(Xcorrected(1:10))/2, 'or', 'LineWidth', 3);
hold off;
grid on;
title('Real value of FFT result', 'FontSize', 18);
xlabel('frequency', 'FontSize', 18);
ylabel('real', 'FontSize', 18);
subplot(3,1,3);
hold on;
stem(f(1:10), pi*imag(X(1:10))/2, 'ob', 'LineWidth', 3);
plot(f(1:10), pi*imag(Xcorrected(1:10))/2, 'or', 'LineWidth', 3);
hold off;
grid on;
title('Imaginary value of FFT result', 'FontSize', 18);
xlabel('frequency', 'FontSize', 18);
ylabel('imag', 'FontSize', 18);

FFT 和锌补偿 FFT 的结果

前十个非零频率的相位值：

>> X_angle(1:10)*2/pi
ans =
-1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000