FFT频谱相频特性仿真结果的问题

guodongliang

用FFT进行相频特性研究，如果横坐标都是模拟频率（Hz），matlab程序如下：
clc
clear
%---------------------- non-uniform sampling
N=1000;
n = 1:N;
f0=100; f1=200; f2=400;% the frequecy of signal
fs=1000; % the sampling frequency
t=n/fs;  
x=sin(2*pi*f0*t)+sin(2*pi*f1*t+pi/2)+sin(2*pi*f2*t+pi/4); % signal
Xk = fft(x);
F=Xk(1:N/2+1); f=fs*(0:N/2)/N;
subplot(311);
plot(t,x);
title('sampled signal')
xlabel('Time / s')
subplot(312);
plot(f,abs(F))
xlabel('Frequency / Hz')
ylabel('DFT magnitude')
subplot(313);
ph = angle(F)/pi*180;%degree
plot(f,ph);
xlabel('Frequency / Hz')
ylabel('DFT Phase / Degree')
% find the phase
[num,val]=sort(abs(F),'descend');
k=4;
   num(1:k)
   val(1:k)
   ph(val(1:k))
结果如图所示：
问题是：相频的结果（如同噪声）对吗？
最后一段代码是查出幅频图上最大4个幅频对应的频点，然后查对应的相位：

ans =
  500.0000  500.0000  500.0000    0.0000

ans =
   201   401   101   400

ans =
   72.0000   99.0000  -54.0000   90.3364

结果与所设定的初相（0、pi/2 、pi/4）不吻合。
请高手指点。

guodongliang

怎么没人回复呢

songzy41

guodongliang 发表于 2011-5-23 15:32

                               
登录/注册后可看大图

用FFT进行相频特性研究，如果横坐标都是模拟频率（Hz），matlab程序如下：
clc
clear

把程序做了一些修改，程序下面列着。主要修改内容为：
1，时间变量t的设置不正确，改为t=(n-1)/fs;，即初始时刻t=0，这样才有可能求出初始相角（t=0时刻的相角）；
2，从angle得到的相角是针对余弦信号的，与正弦信号差90度；
3，因为3个主要频率对应的谱线分别为101，201，401条谱线相接从这几条谱线中读入计算，其答案为（其中直接计算了信号的幅值）
k=101  Am=1.00000   Freq=100.00000   Angle=0.00000
k=201  Am=1.00000   Freq=200.00000   Angle=90.00000
k=401  Am=1.00000   Freq=400.00000   Angle=45.00000
程序为：
clc
clear
%---------------------- non-uniform sampling
N=1000;
n = 1:N;
f0=100; f1=200; f2=400;% the frequecy of signal
fs=1000; % the sampling frequency
t=(n-1)/fs;  
x=sin(2*pi*f0*t)+sin(2*pi*f1*t+pi/2)+sin(2*pi*f2*t+pi/4); % signal
Xk = fft(x);
F=Xk(1:N/2+1); f=fs*(0:N/2)/N;
subplot(311);
plot(t,x);
title('sampled signal')
xlabel('Time / s')
subplot(312);
plot(f,abs(F))
xlabel('Frequency / Hz')
ylabel('DFT magnitude')
subplot(313);
ph = angle(F)/pi*180+90;%degree
plot(f,ph);
xlabel('Frequency / Hz')
ylabel('DFT Phase / Degree')
fprintf('k=101  Am=%6.5f   Freq=%6.5f   Angle=%6.5f\n',abs(F(101))*2/N,f(101),ph(101));
fprintf('k=201  Am=%6.5f   Freq=%6.5f   Angle=%6.5f\n',abs(F(201))*2/N,f(201),ph(201));
fprintf('k=401  Am=%6.5f   Freq=%6.5f   Angle=%6.5f\n',abs(F(401))*2/N,f(401),ph(401));

guodongliang

感谢大师！
还有一个疑问：
这个模拟信号只有3个频率分量，并无噪声，幅度频谱是只有3个点有频率的，也几乎无频谱泄露，为何相频曲线在所有频点处都有相位呢？ 是计算误差？

songzy41

guodongliang 发表于 2011-5-25 11:41

                               
登录/注册后可看大图

感谢大师！
还有一个疑问：
这个模拟信号只有3个频率分量，并无噪声，幅度频谱是只有3个点有频率的，也几 ...

这是由于计算误差造成的。首先计算机数字的表示都是用有限字长表示的，例如pi，我们知着它是一个无限长的无理数，但在计算机中只能用有限长度表示。pi=3.1415926...，而在计算机中只能表示为3.14159，后面的0.0000026...就是有限字长的截断误差了。由pi计算出的三角函数，例如sqrt(2)/2，sqrt(3)/2等等原本也应为无理数，现在也按有限字长表示了，就是数的表示中就有误差，用这些带有误差的数参与运算，使误差累计。当然这些误差很小，信号的幅值计算中都看不出来。但在相位计算出就不一样了，例如FFT后X(i)=1.2e-9+j2.4e-10，它的幅值(模值)就很小，但相角angle(1.2e-9+j*2.4e-10)=0.1974就不小了，所以相频曲线在-pi至pi之间振荡是正常的。

guodongliang

感谢大师！
以后多多指教！

wangyouyi

学习下！