回复 15楼 zhwang554 的帖子
本帖最后由 wdhd 于 2016-6-3 10:55 编辑你有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 编辑 ]
回复 16楼 zhwang554 的帖子
非常感谢！学习中。。。。。[ 本帖最后由 yycc2006 于 2009-4-21 16:10 编辑 ]
回复 16楼 zhwang554 的帖子
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 编辑 ] :@$ :@$ :@P
[ 本帖最后由 yycc2006 于 2009-4-21 15:00 编辑 ]
回复 18楼 yycc2006的帖子
你将程序和一组数据写出来,即如N=240,你什么得到X相位=30.0001,y相位=29.9805
我看看.
[ 本帖最后由 zhwang554 于 2009-4-21 15:19 编辑 ]
回复 20楼 zhwang554 的帖子
close all;clc;clear all;fs=8000;
f1=600;
A=16384;
ph1=30;
load x1.mat
load y1.mat
%正弦信号相位的测量
N=;%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的输入数据
=fcn_apFFT_analysis(y1,N(js));
y2=y(start:stop);
=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 编辑 ]
回复 20楼 zhwang554 的帖子
function= fcn_apFFT_analysis(y,N)win =hanning(N)';
winn=conv(win,win);%apFFT需要卷积窗
win2=winn/sum(winn);%窗归一化
y22=y.*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的相位普
a=a2;
p=p2;
TOAST.exe发不上来，我将数据发上来吧。
我测得是第1200样点的相位
[ 本帖最后由 yycc2006 于 2009-4-21 15:50 编辑 ]
回复 22楼 yycc2006的帖子
你的程序编得顶好.因为我没有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=;%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)+;
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 编辑 ]
回复 23楼 zhwang554 的帖子
读您的帖子很有收获！我总结您的结论如下：!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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 编辑 ]
如何判定噪声的影响
要判定是否是系统噪声的影响,你再取一组数据,如果误差曲线不一样了,那就是噪声的影响,如果误差曲线完全相同,那就不是噪声的影响,再找原因.
你的x和y曲线儿乎完全一样,下图中,兰色为x,红色为y,绿色为x-y,你要求的测相位精度为多高?
你在24楼的结论笫3点3、为了减小噪声在apfft程序中必须加hanning卷积窗。应为3、为了减小泄漏在apfft程序中必须加hanning卷积窗。
[ 本帖最后由 zhwang554 于 2009-4-21 18:44 编辑 ]
回复 25楼 zhwang554 的帖子
相位的精度不需要很高，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 编辑 ]
x和y信号的谱分析比较
对N=240的归1后的x和y作apfft谱分析如下图1图1x和y信号的谱分析比较
可见,x信号中产生许多幅度-50db谐波,这是取样,A/D变换,数据转换造成的,信噪比-100db
y信号中信噪比-50db,将全部谐波淹没
系统的信噪比50db
x信号的相位谱很干净,各谐波相位值各异. y信号的相位谱受噪声干扰.
[ 本帖最后由 zhwang554 于 2009-4-22 10:17 编辑 ]
回复 27楼 zhwang554 的帖子
谢谢您！学习中。。。。y数据和x数据测量程序
下附程序是测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-0019.8828e-001 6.0005e+002
512 3.0006e+0015.6848e-0039.9935e-001 6.0001e+002
256 2.9956e+001 -4.4074e-0029.9608e-001 6.0001e+002
128 3.0396e+0013.9584e-0019.8667e-001 6.0085e+002
240 2.9981e+001 -1.9458e-0029.9465e-001 5.9999e+002 480 2.9855e+001 -1.4516e-0019.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)+;%构成长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=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频率
10243.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+002240 3.0000e+001 1.4475e-004 1.0013e+000 6.0000e+002120 3.0000e+001 1.4472e-004 1.0013e+000 6.0000e+002
结果儿乎全部相同, 表明x是一个fs=8000的近理想的600赫,初相30度余弦函数数据.
[ 本帖最后由 zhwang554 于 2009-4-22 18:56 编辑 ]