请教全相位谱分析问题

w89986581

看过有些介绍全相位的文章，重构输入信号以后，为什么得到的谱峰相位和理论值相差很大！谁能指点一二，谢谢！


tt = 0:1/2000:6;
yy = cos(2*pi*100.4*tt+pi/3)+cos(2*pi*150.6*tt+pi/2)+.1*rand(size(tt));
NFFT = 4096;
yy1 = yy(1:NFFT*2-1);
yy1 = yy1(:);

vecter = [NFFT:-1:1,2:NFFT];
vecter = vecter/NFFT;
for ii = 1:NFFT-1,
    yy2(ii) = vecter(ii)*yy1(ii) + vecter(NFFT+ii)*yy1(NFFT+ii);
end
yy2(NFFT) = vecter(NFFT)*yy1(NFFT);

yy2_fft = fft(yy2,NFFT);
yy2_phase = phase(yy2_fft)*180/pi;
yy2_phase = mod(yy2_phase,360);

结果
yy2_phase(206) = 103
yy2_phase(310) = 64.2
:'( :'( :'( :'(

yangzj

呵呵,这个方法我曾经试过,也曾从理论上去想过,最终还是没能弄明白

yangzj

不知道你是按哪篇文章做的

w89986581

一种减少泄漏的新型谱估计方法
离散频谱的全相位校正法

应该是忽略了什么

yangzj

呵呵，我也很想知道这篇文章是怎么做出来的，你做出来的话能不能交流下。我对这篇文章可以说算是比较熟，当时按作都的方法没能做出来

w89986581

这文章引用率很高，为什么做不出来呢？郁闷~

zhwang554

你的程序基本正确,但有几点错误:
1  加的窗是vecter =[1:NFFT,NFFT-1:-1:1],即三角窗,不是vecter = [NFFT:-1:1,2:NFFT];
2  组成的NTTF阶新序列是yy2=[yy2(NFFT) yy2];而原程序是yy2=[yy2 yy2(NFFT)];应将中间点放在第1项
3 原程序tt从0开始,实测的是(2*NFFT-1)序列中间点的相位,即tt=4095点的相位,新程序给出了该点的理论值,与测量值一致
4  若tt=-NFFT+1:NFFT-1,则测的是中间点tt=0的相位,即初相位,  你改变100.4 到100.55 等,测的初相位不变.
4 先在无噪下测
5 原程序是加三角窗(矩形窗与矩形窗尽卷识),两频率间泄漏影响精度.若用双窗(hanning窗与hanning窗卷识),无噪时精度高.

clear;clf;
tt = 0:1/2000:6;
yy = cos(2*pi*100.4*tt+pi/3)+cos(2*pi*150.6*tt+pi/2);
NFFT = 4096;
yy1 = yy(1:NFFT*2-1);
yy1 = yy1(:);
%vecter =  [1:NFFT,NFFT-1:-1:1];
%vecter = vecter/NFFT;
vecter =  conv(hanning(NFFT)',hanning(NFFT)');
vecter = vecter/NFFT/max(vecter);
for ii = 1:NFFT-1,
    yy2(ii) = vecter(ii)*yy1(ii) + vecter(NFFT+ii)*yy1(NFFT+ii);
end
yy2(NFFT) = vecter(NFFT)*yy1(NFFT);
yy2=[yy2(NFFT) yy2];
yy2_fft = fft(yy2,NFFT);
yy2_phase = phase(yy2_fft)*180/pi;
yy2_phase = mod(yy2_phase,360);
disp('中间点的测量值程序')
yy2_phase(206)
yy2_phase(310)
disp('中间点的理论值')
tt=4095;
mod((2*pi*100.4*tt/2000+pi/3)*180/pi,360)
mod((2*pi*150.6*tt/2000+pi/2)*180/pi,360)

yangzj

原帖由 zhwang554 于 2007-5-6 22:59 发表
你的程序基本正确,但有几点错误:
1  加的窗是vecter =[1:NFFT,NFFT-1:-1:1],即三角窗,不是vecter = [NFFT:-1:1,2:NFFT];
2  组成的NTTF阶新序列是yy2=[yy2(NFFT) yy2];而原程序是yy2=[yy2 yy2(NFFT)];应将中间 ...

呵呵,请教一个问题:
你这里说得到的是中间点的相位,那如果想得到起始点的相位怎么办呢?

zhwang554

可参看"如何准确确定信号中强线谱的相位？"
的讨论,我给w89986581贴子中回答了这个问题.

zhwang554

　　<数字信号全相位谱分析和滤波技术>,电子工业出版社,2009.2
　　一书己出版, 书名就来自本贴题目<请教全相位谱分析问题>,
　　楼主很敏锐, 2007年4月,为了直接用fft谱计算得到的相位, 竟能找到全相位fft,当时全相位fft测相位的论文还没有发表,只是论文<离散频谱的全相位校正法>中提了一下全相位fft中相位不必校正,他竟知全相位fft可直接测相位,找不到代码就自已编, ,成功不了,郁闷,想自已一定忽略了什么,网上发贴子请教,有了回答,”Thanks a million!!!”. 他提出的”全相位谱分析”简称, 以后常用.在google上搜索”全相位”,第一个就是<请教全相位谱分析问题,信号处理方法,振动论坛>.
　　[ 本帖最后由 zhwang554 于 2009-3-17 06:10 编辑 ]

yycc2006

我正在学习<数字信号全相位谱分析和滤波技术>一书，有一个问题向大家请教。
附录1的程序如下：
%FFT和apFFT比较程序，书附录1
close all;clc;clear all;
N=128;
t=-N+1:1:N-1;
y= cos(2*pi*t*(20.4)/N)+0.001*cos(2*pi*t*(28.2)/N+pi);%正弦初相为0和pi
y1 = y(N:2*N-1);%取后N个数据为FFT的输入数据
win =  hanning(N)';
win1 = win/sum(win);%窗归一化
y11= y1.*win1;
y11_fft = fft(y11,N);
a1 = abs(y11_fft);
p1= mod(phase(y11_fft)*180/pi,360);%不用加，没有相位误差
y2=y(1:2*N-1);%取2N-1个数据为apFFT的输入数据
winn=conv(win,win);%apFFT需要卷积窗
win2=winn/sum(winn);%窗归一化
y22=y2.*win2;
y222=y22(N:end)+[0 y22(1:N-1)];%构成长N的FFT输入数据
y2_fft=fft(y222,N);
a2=abs(y2_fft);%apFFT的振幅普
p2=phase(y2_fft)*180/pi;   
p2=mod(p2,360);%apFFT的相位普%
tt=0:N-1;
subplot(211);plot(tt,10*log10(a1),'b.-',tt,10*log10(a2),'r.-');
title('amplitude spectrum');
ylim([-150,30]);xlim([0,N/2]);
legend('fft','apfft');
xlabel('f');ylabel('db');
grid
subplot(212);plot(tt,p1,'bo-',tt,p2,'r.-');
title('phase spectrum');
ylim([0,400]);xlim([0,N/2]);
legend('fft','apfft');
xlabel('f');ylabel('度');
grid

输入信号正弦波的初相位是0和pi，计算apFFT的相位谱p2=mod(phase(y2_fft)*180/pi,360); 正常来说结果应该是0度和180度，但实际结果为 appendix1.fig (17.95 KB, 下载次数: 54)

2009-4-17 14:19 上传

点击文件名下载附件

file:///e:/apFFT/appendix1.jpg
         0
  360.0000
    0.0000
         0
    0.0000
    0.0000
  360.0000
    0.0000
  360.0000
  360.0000
.........
  从而使得相位谱图不为一直线。我已经做过模360度的运算了。请问是什么原因？如何解决？
但如果把初相改成pi/4或其他非0的相位值，结果都正确。只有初相为0时出现上述结果。

[ 本帖最后由 yycc2006 于 2009-4-17 14:29 编辑 ]

zhwang554

                               图1  二个频率成份的fft和apfft对数振幅谱和相位谙
你得到的曲线如上图1.
在第二个频率28.2附近有5个点p2相位是180度左右,没有错.你在11楼取的点少了,再多取几个180度就出来了.
在第一个频率20.4附近p2相位是0度或360度(就是零度)左右,也没有错
在matlab中用长数据读20.4附近到28.2附近的p2值为(在11楼是用短数据读,359.9999都读为360了)
p2(15:35) =
  Columns 1 through 4
    1.503391790030778e-010    3.599999999999483e+002    3.599999999999901e+002    3.323916266600944e-012
  Columns 5 through 8
    3.599999999999996e+002                         0                        1.228441338581712e-016    3.600000000000000e+002
  Columns 9 through 12
    7.820001033112934e-016    3.599999999999998e+002    3.599999999999959e+002    3.599999999999984e+002
  Columns 13 through 16
    8.763350463573120e-011    1.800000000000017e+002    1.800000000000001e+002    1.800000000000000e+002
  Columns 17 through 20
    1.800000000000179e+002    1.799999999920851e+002    3.599999999888563e+002    3.599999999865641e+002
其中不少是359.99999999,不到360,所以模360度的运算以后不在0度处,而在359.99999999度处.
      产生的原因是这儿有二个频率成份,笫2个频率成份28.2虽小,但它在笫一个频率成份附近的泄漏使其相位值偏离0度,为359.9999999,从而曲线在0和360跳变.你可用-200度到+200度的Y轴代现在的0到400度Y轴,这时0到360度的跳变就不存在了(但这时正负180度有一跳变)
      如果你将28.2的振幅值置零,可得下图2的一个频率成份的水平相位特性.所以全相位水平相位特性严格说来只有单频率成份能成立
      虽然一个频率成分峰值附近的相位值差不多,但只有振幅谱峰值处的相位值才是这个频率成分的相位值,特别是有噪时更应用峰值处的相位值

                             图2  一个频率成份的fft和apfft对数振幅谱和相位谙

[ 本帖最后由 zhwang554 于 2009-4-17 20:24 编辑 ]

zhwang554

    喔,我明白了.你在学习<数字信号全相位谱分析和滤波技术>一书时,重复了书附录1的程序 %FFT和apFFT比较程序.
    发现当信号初相0和pi时
y= cos(2*pi*t*(20.4)/N)+0.001*cos(2*pi*t*(28.2)/N+pi);
    Apfft的相位谱在0和360度之间跳变,如12楼中图1所示,不能清晰反映apfft的水平相位特性,引起初学者迷惑.

    应将初相改为pi/2和pi
y= cos(2*pi*t*(20.4)/N+pi/2)+0.001*cos(2*pi*t*(28.2)/N+pi)
     这样apfft的相位谱如下图3,充分反映apfft的水平相位特性.

           图3  二个频率成份的fft和apfft对数振幅谱和相位谱
     
     这是书中考虑不周之处.以后有机会再引用附录1的FFT和apFFT比较程序,一定把它改过来,便於初学者理解.谢谢你的指点.很高兴你学习<数字信号全相位谱分析和滤波技术>一书,有问题请提出.




[ 本帖最后由 zhwang554 于 2009-4-20 08:08 编辑 ]

yycc2006

非常感谢您的解答！讲解得很透彻！明白啦！:lol
我把初相改为pi/2和pi时，得到了预期的结果。:victory:

另外，我正尝试将全相位FFT应用于经系统传输引起失真的单频正弦信号的频谱分析，目前还有一些问题。
等仿真结果出来再向您请教！

[ 本帖最后由 yycc2006 于 2009-4-20 10:09 编辑 ]

yycc2006

FFT的长度确定问题
我将全相位FFT应用于经系统传输引起失真的单频正弦信号的频谱分析，获得了初步结果！
但是在确定FFT长度时不知道该采用什么准则？
1、需解决的问题：考察一个黑箱系统对0-4KHz频率的响应,系统的输入要求采样率fs=8000Hz
2、解决的过程：
（1）以50Hz步进选择信号频率：50、100、150.....
（2）构造余弦函数x=cos(2pi*t*f /fs+30*pi/180);
（3）经系统传输后得到输出y；
（4）对y进行频谱分析，分别采用FFT和apFFT方法
以f=800Hz为例，apFFT的幅度谱为：（N=2000）

                                    图1 apFFT分析所得幅度谱
相位谱为：

                                 图2  apFFT分析所得相位谱
从图2上看相位谱很乱，没有表现出“水平直线特性”，但是以谱峰值点为中心展开后发现，在谱峰值点201位置附近有3点近似水平特性，见图3.

                      图3 apFFT分析所得相位谱（展开）

根据以上仿真结果，我认为：
（1）输出信号y的频率F=800Hz，没有产生新的频率，失真是线性的；
（2）输出信号y的相位发生了偏移；
（3）输出信号y的幅度有失真（这可以直观地从时域波形上看出）。
现在的问题是：相位偏移失真有多大？
在进行频谱分析时发现，不同长度的FFT，所测得的相位是不同的。如果N选择不合适（如N=600)，apFFT所测结果将不如FFT的精确。如图4.                                    

                           图4  FFT和apFFT在不同分析长度下测得的初相位
请问：
（1）不同的N测出了不同的相位，我应该选择哪个相位值作为信号y的初相位？
（2）在进行apFFT谱分析时，该如何来确定FFT的长度N？
（3）当采样率fs不是f的整数倍时（如f=300Hz）,N如何确定？
（4）当f值增大（如f=3500Hz），一个周期内的采样点数将减少。进行apFFT分析时要不要进行时域插值？精度够吗？
或者说，进行apFFT分析时，需要有多高的采样率？



（有谁知道附件怎么删除啊？不小心上载错了，呵呵）

[ 本帖最后由 yycc2006 于 2009-4-20 12:50 编辑 ]