w89986581 发表于 2007-4-28 11:31

请教全相位谱分析问题

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


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 = ;
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 发表于 2007-4-28 13:07

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

yangzj 发表于 2007-4-28 13:08

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

w89986581 发表于 2007-4-28 14:26

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

应该是忽略了什么

yangzj 发表于 2007-4-28 14:35

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

w89986581 发表于 2007-4-28 18:48

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

zhwang554 发表于 2007-5-6 22:59

全相位FFT测相位

你的程序基本正确,但有几点错误:
1加的窗是vecter =,即三角窗,不是vecter = ;
2组成的NTTF阶新序列是yy2=;而原程序是yy2=;应将中间点放在第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 =;
%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_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 发表于 2007-5-7 10:44

本帖最后由 wdhd 于 2016-6-3 10:54 编辑

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

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

zhwang554 发表于 2007-5-7 12:29

回复 #8 yangzj 的帖子

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

zhwang554 发表于 2009-3-17 05:38

全相位谱分析

本帖最后由 wdhd 于 2016-6-3 10:54 编辑

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

yycc2006 发表于 2009-4-17 14:13

回复 10楼 zhwang554 的帖子

我正在学习<数字信号全相位谱分析和滤波技术>一书,有一个问题向大家请教。
附录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)+;%构成长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();
legend('fft','apfft');
xlabel('f');ylabel('db');
grid
subplot(212);plot(tt,p1,'bo-',tt,p2,'r.-');
title('phase spectrum');
ylim();xlim();
legend('fft','apfft');
xlabel('f');ylabel('度');
grid

输入信号正弦波的初相位是0和pi,计算apFFT的相位谱p2=mod(phase(y2_fft)*180/pi,360); 正常来说结果应该是0度和180度,但实际结果为
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 发表于 2009-4-17 16:27

回复 11楼 yycc2006 的帖子


                               图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 发表于 2009-4-20 06:18

再回复 11楼 yycc2006 的帖子

    喔,我明白了.你在学习<数字信号全相位谱分析和滤波技术>一书时,重复了书附录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 发表于 2009-4-20 10:06

回复 13楼 zhwang554 的帖子

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

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

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

yycc2006 发表于 2009-4-20 12:15

回复 13楼 zhwang554 的帖子

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分析所得幅度谱
相位谱为:

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

                      图3 apFFT分析所得相位谱(展开)

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

                           图4FFT和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 编辑 ]
页: [1] 2 3
查看完整版本: 请教全相位谱分析问题