如何准确确定信号中强线谱的相位？

zhwang554

1。你执行n次程序是为了什么，如果是测噪声下多次测量相位的结果，你的程序不对，每执行一次y22不会变。

2 还有如果按照全相位的定义，
for ii = 1:NFFT-1,
    yy2(ii) = vecter(ii)*yy1(ii) + vecter(NFFT+ii)*yy1(NFFT+ii);
end
是否应该改为：
for ii = 1:NFFT-1,
    yy2(ii) = i*vecter(ii)*yy1(ii) + (NFFT-i)*vecter(NFFT+ii)*yy1(NFFT+ii);
end

vecter中已经包含  i  和 (NFFT-i)，不必再加权。

3。最后能否解释一下vecter = vecter/NFFT/max(vecter);的意思？

vecter是卷积窗函数，上式是卷积窗函数归一。

yifeng2030

1.我循环的本意是要测噪声下多次测量相位的结果，当然这个程序不完整，我发现加一个100次的循环后，yy2数组长度确实是1*4196；你若有空的话，请用matlab试试
clear all;
T=100;
while T>0
   T=T-1;
tt = -4095/2000:1/2000:4095/2000;
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)
end



2。按照全相位的定义，可否改为：
vecter = vecter/max(vecter);%归一化
for ii = 1:NFFT-1,
    yy2(ii) = i*vecter(ii)*yy1(ii) + (NFFT-i)*vecter(NFFT+ii)*yy1(NFFT+ii);
end
yy2(NFFT) = NFFT*vecter(NFFT)*yy1(NFFT);
yy2=[yy2(NFFT) yy2]/NFFT;
这样的话是否更符合全相位的定义？
还有，你有关于频谱插值校正法的材料吗？有的话，可否发给我?yifeng2030@126.com，不胜感激。
呵呵，有点罗索。。。。。:@)

couldnt2000

:victory: :victory: :victory:

zhwang554

你的问题有道理,主要是原程序中变量yy1(ii)和yy2(NFFT)在循环语句中不是一个数,而是不断累积,虽然FFT时只取NFFT个不影响计算结果,但yy2还是增长了个数,
多次测量用循环语句可用下程序,避克出现变量yy1(ii)及yy2(NFFT),更快,更简单

clear all;close all
NFFT = 4096;
T=10;
while T>0
   T=T-1;
tt = -(NFFT-1)/2000:1/2000:(NFFT-1)/2000;
yy =2* cos(2*pi*100.4*tt+pi/3);
yy1 = yy(1:NFFT*2-1);
vecter =  conv(hanning(NFFT)',hanning(NFFT)');
vecter = vecter/sum(vecter);
yy1a= yy1.*vecter;
yy2=yy1a(NFFT:end)+[0 yy1a(1:NFFT-1)];
yy2_fft = fft(yy2,NFFT);
yy2_phase = phase(yy2_fft)*180/pi;
yy2_phase = mod(yy2_phase,360);
disp('初相位测量值')
yy2_phase(round(100.4*NFFT/2000)+1)
end

2。按照全相位的定义，可否改为：
vecter = vecter/max(vecter);%归一化
for ii = 1:NFFT-1,
    yy2(ii) = i*vecter(ii)*yy1(ii) + (NFFT-i)*vecter(NFFT+ii)*yy1(NFFT+ii);
end
yy2(NFFT) = NFFT*vecter(NFFT)*yy1(NFFT);
yy2=[yy2(NFFT) yy2]/NFFT;
这样的话是否更符合全相位的定义？

上公式中的 i 是 ii ?
在原程序中,vecter是卷积窗,当二个矩阵窗作卷积时,生三角窗,
三角窗vecter(ii)等於ii,vecter(NFFT+ii)等於(NFFT-ii),所以上面的式中就不必再乘 i 和 (NFFT-i)了.
归一指输入信号幅值为1,最大频谱振幅值也是1,矩形窗用vecter = vecter/max(vecter)/NFFT或任意窗用vecter = vecter/sum(vecter)

全相位信号形成
1 . 取2N-1个信号                                                            yy1=yy(1:2*N-1)
2.  用两个长N的对称窗产生长2N-1归一卷积窗                      vecter
3.  长2N-1信号乘长2N-1卷积窗                                          yy1a=yy1.*vercter
4.  长2N-1 的yy1a移位相加形式长N的全相位输入数据 yy2  
       即yy1a的前N-1个数据补0形成N个数数据                                                [0   yy1a(1:N-1)]
       和yy1a的后N个数琚                                                         yy1a(N;end)
       两者相加形成长N的全相位输入数据                             yy2=yy1a(N:end)+[0  yy1a(1:N-1)]

[ 本帖最后由 zhwang554 于 2007-5-26 14:25 编辑 ]

yifeng2030

1 对,i写错了,就是ii
呵呵，终于看明白了，但感觉和fft的精度差不多，有没有什么方法可以再提高精度？

[ 本帖最后由 yifeng2030 于 2007-5-26 10:29 编辑 ]

zhwang554

全相位FFT的泄漏是原FFT泄漏的平方,泄漏分贝减少一倍,如FFT泄漏-20db,全相位
FFT泄漏为-40db.你用apFFT的程序将apFFT的振幅谱和原FFT的振幅谱画出来,一比就知道了.

在无噪时,全相位FFT测相拉明显好於FFT,你测时发现了幺?

但FFT泄漏是很小的,所以全相位FFT在几乎无噪时才明显好於FFT.但在密集频谱时,互相泄漏影响大,
全相位FFT好些.
    有噪声时,理论上全相位FFT只比传统FFT的精度好2/3,实用优势不如无噪,但总比FFT
好些.但在有噪相位测量时,apFFT由於不用校正,一次FFT测相位,精度虽类同,但全相位
FFT测相位比传统FFT始终有实用优势.
所以全相位测相位无噪时好一倍,有噪时略好些,密集谱好多些.

FFT的相位你什么测出来?一次FFT测得出来吗?校正后?

[ 本帖最后由 zhwang554 于 2007-5-26 11:00 编辑 ]

yifeng2030

我做的都是一百次后求平均精度，没有校正，apfft作N点FFT精度比2N-1点FFT的精度差不了多少，能留个QQ吗？我的是95424388，有问题方便向您请教，呵呵:@) :@)

zhwang554

很高兴和你交流,
95224388是什么?

没有校正,作FFT相位什么测出来?FFT测相位不是有误差呜?和频偏有关?如上述程序中60度初相位FFT能一次FFT测出吗?

你说"apfft作N点FFT精度比2N-1点FFT的精度差不了多少",那么2N-1点FFT的计算量就大了一倍.

而且2N-1点FFT的泄漏比apFFT作N点的泄漏还大.你作一下实验可知.

这问题在本论坛"离散频谱校正"中讨论过,有附图,可参看.

[ 本帖最后由 zhwang554 于 2007-5-26 16:51 编辑 ]

eight

原帖由 zhwang554 于 2007-5-26 11:18 发表
95224388是什么,电话号码?什么打?

QQ...........

yifeng2030

请问楼上两位：
我要对有限长的两个正弦信号通过相关求相位差，那我应该是求圆周相关还是线性相关？能否提供按定义求相关的matlab程序？

zhwang554

你在#17提的原程序的问题,再仔细分析一下,只要稍作改动,原程序可以循环使用
1  将 yy2(NFFT) = vecter(NFFT)*yy1(NFFT)  改为 yy2a = vecter(NFFT)*yy1(NFFT);改为一定量
2  yy2变量多次使用,在for-end中使用一次,在yy2=[yy2(NFFT) yy2]中又使用一次;  改为另一变量yy2b=[yy2a yy2];这是主要问题.一次使用
没有问题,循环使用就累积了.
3  由於窗是hanning(NFFT)的卷积,归一用sum(vecter),这样振幅为2的余弦信号在无频偏时振幅最大峰值为1.
这个程序原由w89986581辛苦钻研后编成,我在保留原貌下改了一下,你循环使用发现一些问题,改后又进一步.
修改后的程序

close all;
clear all;
T=4;
k=1;
while T>0
NFFT = 4096;
tt=(-NFFT+1)/2000:1/2000:(NFFT-1)/2000;
yy =2* cos(2*pi*100.4*tt+pi/3)+0.0*rand(size(tt));
yy1 = yy(1:NFFT*2-1);
yy1 = yy1(:);
vecter =  conv(hanning(NFFT)',hanning(NFFT)');
vecter = vecter/sum(vecter);
for ii = 1:NFFT-1,
   yy2(ii) = vecter(ii)*yy1(ii) + vecter(NFFT+ii)*yy1(NFFT+ii);
end
yy2a = vecter(NFFT)*yy1(NFFT);
yy2b=[yy2a yy2];
yy2_fft = fft(yy2b,NFFT);
yy2_phase = phase(yy2_fft)*180/pi;
yy2_phase = mod(yy2_phase,360);
phase_1(k)=yy2_phase(round(100.4*NFFT/2000)+1)
T=T-1;
k=k+1;
end

[ 本帖最后由 zhwang554 于 2007-5-27 13:08 编辑 ]

yifeng2030

在原程序中,vecter是卷积窗,当二个矩阵窗作卷积时,生三角窗,
三角窗vecter(ii)等於ii,vecter(NFFT+ii)等於(NFFT-ii),所以上面的式中就不必再乘 i 和 (NFFT-i)了.
归一指输入信号幅值为1,最大频谱振幅值也是1,矩形窗用vecter = vecter/max(vecter)/NFFT或任意窗用vecter = vecter/sum(vecter)



那么如果w=conv(hanning(NFFT)',hanning(NFFT)'),卷积窗就不是三角窗了
for ii = 1:NFFT-1,
   yy2(ii) = vecter(ii)*yy1(ii) + vecter(NFFT+ii)*yy1(NFFT+ii);
end
要改为：
for ii = 1:NFFT-1,
    yy2(ii) = i*vecter(ii)*yy1(ii) + (NFFT-i)*vecter(NFFT+ii)*yy1(NFFT+ii);
end
么？

zhwang554

那么如果w=conv(hanning(NFFT)',hanning(NFFT)'),卷积窗就不是三角窗了
for ii = 1:NFFT-1,
   yy2(ii) = vecter(ii)*yy1(ii) + vecter(NFFT+ii)*yy1(NFFT+ii);
end
要改为：
for ii = 1:NFFT-1,
    yy2(ii) = i*vecter(ii)*yy1(ii) + (NFFT-i)*vecter(NFFT+ii)*yy1(NFFT+ii);
end
么？

不要改,上面的公式是通式,事实上全相位加窗有三种情况(从全相位FIR滤波器中引入,对称窗指hann等常用窗)
无窗  二个矩形窗                           卷积
单窗  一个矩形窗,一个对称窗           卷积
双窗  二个对称窗                           卷积

在频谱分析中,对单指数信号,全相位FFT无窗的泄漏   是 FFT无窗泄漏的平方
                                      全相位FFT双窗的泄漏   是 FFT加窗泄漏的平方
所以常用无窗,双窗
相位测量和窗无关,但旁频泄漏影响精度,全相位双窗泄漏最小(FFT加窗后泄漏比FFT不加窗小了,全相位双窗泄漏比加窗FFT泄漏分贝更小一半),所以
                         测相位必用双窗测,因为就单余弦波,其镜像复指数谱线的泄漏就影响原指数谱线的相位测量

它们统一用上面的式子即可,
即 无窗        vecter=conv(ones(1,NFFT),ones(1,NFFT))
或 双窗        vecter=conv(hanning(NFFT)',hanning(NFFT)')
统一用
for ii = 1:NFFT-1,
   yy2(ii) = vecter(ii)*yy1(ii) + vecter(NFFT+ii)*yy1(NFFT+ii);
end

如果输入数据和卷积窗先乘好了
yy1a=yy1.*vecter
则循环更简单,更易理解
for ii = 1:NFFT-1,
   yy2(ii) = yy1a(ii) + yy1a(NFFT+ii);
end
当时为了尽量保留原程序语句,没有先乘

[ 本帖最后由 zhwang554 于 2007-5-28 18:18 编辑 ]

yifeng2030

呵呵，多谢了

zhwang554

利用w89986581的测相位程序(#26)可以作许多实验:

1 将信号相位改写为
    yy =2* cos(2*pi*100.4*tt+60*pi/180);
这时输出是60度
将初相位改写为
    yy =2* cos(2*pi*100.4*tt+60.1111111111111*pi/180);
这时输出是60.1111111111108度
这说明无噪全相位FFT测相位是很准的,达10^-11,(校正后FFT可这60.111111,达10^-6,差不多一半,和apFFT泄漏分贝比FFT小一倍一致)

你可以加一个或多个信号
    yy =2* cos(2*pi*100.4*tt+60.11111111111111*pi/180)+2*cos(2*pi*150.6*tt+90.1111111111111*pi/180);
这时输出是60.1111111111138度
               90.1111111111078度
相位测量仍很准,互相影响不大,因为apFFT泄漏小,互相影响小

2 信号中的频率你什么变, 如100.4 改为100.3  100.2
这时输出是60.1111111111138度,  60.1111111111144度,  60.1111111111159度
相位测量不变,
我们是测谱线峰值处的相位,不管频率偏离多少,在这一点测相位仍对,这个性质有点不好理解,在离散点可测频率连续变化时的相位,
但实测是如此,理论公式也证明如此,用全相位信号组成也解释得通,使DFT(离散Fourier变换)特别有用
全相位FFT的相位和频率偏离无关

3 你画出相位特性
subplot(2,1,1);stem(abs(yy2_fft),'.');
subplot(2,1,2);plot(yy2_phase);
(见附图)
相位特性在谱峰值附近是一条水平线,即频率偏离超过频谱线1格以上,相位测量仍正确,
如果是单复指数信号,,apFFT的相位特性是一条水平线(处处等於初相位值),FFT的相位特性是一条斜线(频偏时在峰值谱线处的相位也不等於初相位值)
当然,有噪时还是峰值处测的相位值最好

4你再加噪音,多次测均方差值,均不比FFT差,但优势不大,密集频谱好多些
  要和fft作相位测量比较麻烦,因为fft的相位必用频偏校正,所以还要将频偏较正过来才行.

[ 本帖最后由 zhwang554 于 2007-5-29 07:17 编辑 ]