请教全相位谱分析问题

zhwang554

　　你有x,y两组数据,同时取2N-1个数据,分别对x和y作apfft,比较振幅和相位.
　　apfft测的是2N-1个数据中间点的相位,对x数据这点相位不会是初相位30度,它和t等有关,但这没有关系,同时测出x和y的相位,相位差就知道了,就知道x和y有没有相移.
　　这和用双线示波器观察x和y二个波形一样,但示波器不能精确读出振幅,相位和频率的数值.
　　用FFT也可以读出振幅和频率,但相位只能读出相位差,不能正确读出x和y的瞬间相位(除非经过麻频的校正).
　　而apfft可以分别读出x和y的瞬间相位,二者一减相位差也出来了.
　　你的问题可参见下贴子
　　<关于离散信号的幅度和相位估计>
　　作者: grassing 时间: 2009-3-9 16:15 研学论坛
http://bbs.matwav.com/viewthread.php?tid=854656&extra=page%3D3




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

yycc2006

非常感谢！学习中。。。。。

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

yycc2006

zhwang554 您好！
对输入x,输出y两组数据,同时取2N-1个数据,分别对x和y作apfft,比较振幅和相位.分别读出x和y的瞬间相位,二者相减获得相位差。
但是我在仿真中发现，分析的长度N不同，其相位差的值是不同的，而且差别很大。
对于cos(2*pi*t*f1/fs+ph1*pi/180)   其中f1=600Hz       fs=8000
因为80和6的最小公倍数是240，因此N取240的倍数。
以下是仿真结果：
长度N：  240            360         480           720         960       1200
X相位：30.0001   30.0001   30.0001   30.0001   30.0001   30.0001
y相位：29.9805   29.8849   29.8548   29.8424   29.8645   29.9209
相差： -0.0196    -0.1152     -0.1453   -0.1577    -0.1357    -0.0793
其中，N=240时，获得最小的相差
因为我的目的是对某黑箱系统进行辨识，因此需要精确地求出相移。
请教：如何来确定N的值呢？采用试验方法吗？谢谢！

[ 本帖最后由 yycc2006 于 2009-4-21 15:03 编辑 ]

yycc2006

:@$ :@$ :@P

[ 本帖最后由 yycc2006 于 2009-4-21 15:00 编辑 ]

zhwang554

你将程序和一组数据写出来,
即如N=240,你什么得到X相位=30.0001,y相位=29.9805   
我看看.

[ 本帖最后由 zhwang554 于 2009-4-21 15:19 编辑 ]

yycc2006

close all;clc;clear all;
fs=8000;
f1=600;
A=16384;
ph1=30;
load x1.mat
load y1.mat
%正弦信号相位的测量
N=[240 360 480 720 960 1200 ];%6
for js=1:6
                        middle=1200;%x1.mat,y1.mat数据时求第1200点的相位
                        start=middle-N(js)+1;
                        stop=middle+N(js)-1;
                        y1 =x(start:stop);%取2N-1个数据为apFFT的输入数据
                        [a1 p1]=fcn_apFFT_analysis(y1,N(js));                        
                        y2=y(start:stop);
                        [a2 p2]=fcn_apFFT_analysis(y2,N(js));
                        r=round(f1/fs*N(js))+1;%谱峰值所处的谱线
                        pppx(js)=p1(r);              
                        pppy(js)=p2(r);                                                   
                        d_p(js)=p2(r)-p1(r);  
end%js
d_p
pppx
pppy
plot(d_p,'.-');

[ 本帖最后由 yycc2006 于 2009-4-21 16:02 编辑 ]

yycc2006

function  [a,p]= fcn_apFFT_analysis(y,N)

        win =  hanning(N)';
        winn=conv(win,win);%apFFT需要卷积窗
        win2=winn/sum(winn);%窗归一化
        y22=y.*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的相位普
        
        a=a2;
        p=p2;
TOAST.exe发不上来，我将数据发上来吧。
x1.mat (309 Bytes, 下载次数: 40)

2009-4-21 15:46 上传

点击文件名下载附件

y1.mat (3.95 KB, 下载次数: 35)

2009-4-21 15:46 上传

点击文件名下载附件

我测得是第1200样点的相位

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

zhwang554

你的程序编得顶好.
因为我没有y的数据,我下面的程序是只求x的
close all;clc;clear all;
N1=160;
fs=8000;
t=0:1:300*N1-1;%取300个帧
f1=600;
A=16384;
ph1=30;
x= A*cos(2*pi*t*f1/fs+ph1*pi/180);
%正弦信号频率、幅度和相位的测量
N=[240 360 480 720 960 1200 ];%6
for js=1:6
                        middle=10*160+1;%求start1点的相位,从1601点开始
                        start=middle-N(js)+1;
                        stop=middle+N(js)-1;
                        y1 =x(start:stop);%取2N-1个数据为apFFT的输入数据
win=hanning(N(js))';win1=triang(N(js))';
win2=conv(win,win);win1=win/sum(win);
win2=win2/sum(win2);
s2=y1(1:2*N(js)-1);%第1组(2N-1)个数据
      y2=s2.*win2;
      y2a=y2(N(js):end)+[0 y2(1:N(js)-1)];
      Out2=fft(y2a,N(js));
      a2=abs(Out2);
      p2=mod(phase(Out2),2*pi);
     r=round(f1/fs*N(js))+1
     ppp=p2(r)*180/pi         
end
运行结果:
N=240      ppp =29.9999999995375
N=360      ppp =29.9999999999584
N =480     ppp =29.999999999992
N=720      ppp =29.9999999999986
N=960      ppp =29.9999999999991
N=1200    ppp =29.999999999999
和你的都是30.00001还有些差别?正常呜?
你问如何选N, 上面的相位测量精度随N而提高.,这是对的,因为你的信号是余弦,它有一个镜像频率的泄漏影响精度, N越大,镜像频率越运,泄漏影响越小,
另外噪声影响也是N越大越好,
这和你的实验结论不同.
另外apfft中信号不必整数倍取样,N任何值都可以.

你的系统有噪声呜?如果有y中的噪声影响就在0.2度左右
看来你的系统x和y相移很小. 所以必减小噪声. apfft程序中必须加hanning卷积窗

[ 本帖最后由 zhwang554 于 2009-4-21 16:45 编辑 ]

yycc2006

读您的帖子很有收获！我总结您的结论如下：
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1、相位测量精度随N而提高。因为余弦信号有一个镜像频率的泄漏影响精度, N越大,镜像频率越运,泄漏影响越小,
2、减小噪声影响也是N越大越好，
3、为了减小泄漏在apfft程序中必须加hanning卷积窗。
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

我将数据上传在22楼。

“和你的都是30.00001还有些差别?正常呜?”
我想30.00001结果的原因，可能是因为：对于x=cos，x是double类型的，但信道输入输出要求是16位的整型，因此在数据转换时带来了误差。为了使得输入输出的精度保持一致，在对x求apFFT时，我没有用double型的x，而是用int16型的x来求的。这可能就是得到30.00001，与您的结果29.99999有差别的原因。

“这和你的实验结论不同.”
对于接收端数据y，仿真结果显示精度没有随着N的增加而提高。:funk::@L

[ 本帖最后由 yycc2006 于 2009-4-21 19:21 编辑 ]

zhwang554

要判定是否是系统噪声的影响,你再取一组数据,如果误差曲线不一样了,那就是噪声的影响,
如果误差曲线完全相同,那就不是噪声的影响,再找原因.

你的x和y曲线儿乎完全一样,下图中,兰色为x,红色为y,绿色为x-y,

你要求的测相位精度为多高?

你在24楼的结论笫3点

3、为了减小噪声在apfft程序中必须加hanning卷积窗。

应为

3、为了减小泄漏在apfft程序中必须加hanning卷积窗。

[ 本帖最后由 zhwang554 于 2009-4-21 18:44 编辑 ]

yycc2006

相位的精度不需要很高，10e-4就行了。只是因为我做不同长度N的apFFT得到的相位值不同，所以不知道该取哪一个是正确反映接收端数据y的相位值。判别的依据可以用min|Phasey-Phasex|吧？
由于Y数据的相位精度没有随着N的增加而提高，因此我现在是用试验的方法，通过改变N值来计算Y的相位，找出相差最小的那个相位。
由于需要测试的频率分别为50、100、150.。。。4000，对不同的频率f，fs/f的值不同，或为整数倍或不为整数倍，这样就得采用不同的N去试验。工作量不小，而且没能从理论上说明白。

[ 本帖最后由 yycc2006 于 2009-4-21 21:24 编辑 ]

zhwang554

对N=240的归1后的x和y作apfft谱分析如下图1

                         图1  x和y信号的谱分析比较

可见,x信号中产生许多幅度-50db谐波,这是取样,A/D变换,数据转换造成的,信噪比-100db
        y信号中信噪比-50db  ,将全部谐波淹没
系统的信噪比50db
x信号的相位谱很干净,各谐波相位值各异. y信号的相位谱受噪声干扰.

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

yycc2006

谢谢您！学习中。。。。

zhwang554

    下附程序是测y数据的频率,振幅和相位.是书中附录2,fft/apfft校正程序.这个方法可参看书中 笫4.3.4节 FFT/apFFT 综合相位差校正法.
    这样在告之fs(发送端取样时保证数据中点相位是初相位)下,接收端用N任意值都可测y的3个参数.接收端不必有x数据,认为是600赫,初相位30度,振幅1的余弦信号.

结果和下:

     N      y相位        y相位误差       y振幅         y频率

   1024   2.9770e+001 -2.2972e-001  9.8828e-001     6.0005e+002

    512   3.0006e+001  5.6848e-003  9.9935e-001     6.0001e+002

    256   2.9956e+001 -4.4074e-002  9.9608e-001     6.0001e+002

    128   3.0396e+001  3.9584e-001  9.8667e-001     6.0085e+002

    240   2.9981e+001 -1.9458e-002  9.9465e-001     5.9999e+002

    480   2.9855e+001 -1.4516e-001  9.9463e-001     6.0003e+002

close all;clc;clear all;
fs=8000;
f1=600;
A=16384;
ph1=30;
load x1.mat
load y1.mat
N=256 ;
middle=1200;%x1.mat,y1.mat数据时求第1200点的相位
start=middle;
stop=middle+N-1;
y1 =y(start:stop);%取N数据为FFT的输入数据
y1=y1/max(y1);
start=middle-N+1;
stop=middle+N-1;               
y2=y(start:stop);%取2N-1个数据为apFFT的输入数据
y2=y2/max(y2);
win =  hanning(N)';
win1 = win/sum(win);%窗归1
y11= y1.*win1;
y11_fft = fft(y11,N);
a1 = abs(y11_fft);%FFT振幅谱
p1 = mod(phase(y11_fft)*180/pi,360);;%FFT相位谱
win =  hanning(N)';;
winn =  conv(win,win);%apFFT须要卷积窗
win2 = winn/sum(winn);%窗归1
y22= y2.*win2;
y222=y22(N:end)+[0 y22(1:N-1)];%构成长N的apFFT输入数据
y2_fft = fft(y222,N);;
a2 = abs(y2_fft);%apFFT振幅谱
p2=mod( phase(y2_fft)*180/pi,360);%apFFT相位谱
ee=mod((p1-p2)/180/(1-1/N),1);%频率偏离校正值
aa=(a1.^2)./a2*2;%振幅校正值
r=round(f1/fs*N)+1;%谱峰值所处的谱线
r1=floor(f1/fs*N)+1;%计频率值用
pppx=p1(r);              
pppy=p2(r) ;                                                
d_p=p2(r)-30 ;
aaa=aa(r);
df=mod(ee(r),1);%频偏值
fff=(r1-1+df)*fs/N;
%fff1=(r1-2+df)*fs/N%(N=240);
[N,pppy, d_p,aaa,fff]

上程序中当N=240等整数值取样时, 有时计算频率要减1.即用fff式, 有时计算频率要减2.即用fff1式, 还没能统一为一个式子. 这个问题和你一开始测初相位零时发生0和360度跳变一样, 当理论频偏值为0的信号, 用附录2的程序统一公式正确, 但像y实测数据, 频偏可能是0.9999或0.0001, 所以有时要多减1. 要加一个判断语句即可.
现用的y数据除了N=120,240要用fff1, 其它N=360,480等用用fff. 其它N值(不管df是>0.5或<0.5)都用fff.

   上程序将y改为x也可测x的3个参数,结果为(N任何值时频率计算都用fff).

  N      x相位        x相位误差       x振幅         x频率

1024  3.0000e+001   1.4475e-004   1.0014e+000   6.0002e+002

512   3.0000e+001   1.4475e-004   1.0014e+000   6.0000e+002

256   3.0000e+001   1.4475e-004   1.0014e+000   6.0000e+002

360   3.0000e+001   1.4475e-004   1.0014e+000   6.0000e+002

240   3.0000e+001   1.4475e-004   1.0013e+000   6.0000e+002

120   3.0000e+001   1.4472e-004   1.0013e+000   6.0000e+002

结果儿乎全部相同, 表明x是一个fs=8000的近理想的600赫,初相30度余弦函数数据.



[ 本帖最后由 zhwang554 于 2009-4-22 18:56 编辑 ]

yycc2006

谢谢您！努力学习中。。。。