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

w89986581

如何准确确定信号中强线谱的相位？求fft变换谱的相位？跳变怎么处理的？

[ 本帖最后由 w89986581 于 2007-4-26 17:41 编辑 ]

yangzj

你说的跳变指的是什么

w89986581

呵呵,就是直接用fft谱计算得到的相位,我查了以下有介绍全相位校正的,有没有相关代码可以利用？

这里的蓝线对应谱峰,即希望得到的相位频率，红线对应相位值（度）

[ 本帖最后由 w89986581 于 2007-4-26 18:17 编辑 ]

yangzj

要得到准确的相位信息提取出来，首先需要把准确的频率校正出来。具体的找些论文看下，有比值（插值）法，能量重心法，相位差法等，你说的全相位差法你也可以试试看

zhwang554

你可参看回复 #3 w89986581 的帖子.
下面是全相位FFT测初相位的程序
clear;clf;
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)

w89986581

Thanks a million!!!
构造信号时，为什么时间序列是这样构造的呢？对于一个实际信号，怎么确定相位呢，或者相对相位差？

w89986581

非常感谢Professor Wang！

zhwang554

apFFT测相位需要2N-1个样点，测出的2N-1个样点中间点的时刻的相位值。原FFT测相
位需要N个样点，测出的是第1个样点时刻的相位值,但它须校正.
    实用时,我们对一个正弦波连续取样2N-1点，测出的是第N个样点的相位，但若要知道
第一个取样时刻的样点相位，即初相位,须知第N个样点的时间T,从测量值减去相隔T的相
位值,即初相位.
    但实际上,测一个正弦波的相位,我们不知道它什么时候开始的,测时离起始时间多运.
但这没有关系.
    apFFT测的是任一样点时刻的相位,即样点相位。目前流行的相位计测的都是比较相位，
被测信号和一个同一频率参考信号的相位差。
    相位差指两个同一频率的正弦波的相位差,测出同一时刻的二个正弦波的相位,其差值
就是相位差,任何时刻测出的同一频率二个正弦波的相位差都是一样的.所以同一频率
的两个正弦波的相位差物理意义十分清楚。
    apFFT测相位差就是分别测二个信号在任何同一时刻的相位,其差值即相位差,即同时对
二个同频信号分别取2N-1点,用apFFT测出中间样点相位,其差值即相位差.流行的
相位计直接测被测信号和一个同一频率参考信号的相位差,它不能测样点相位,是十分不同
的.
    样点相位测量是十分有用的,如两个电网要并网,需要同频同相同幅,测出同一时刻的
两路的样点相位是首要的.



apFFT测相位需取2N-1个样点,若取样间隔有变化,apFFT仍正确测量;

你原写的程序中,N不必取4096,取大的好处是有噪声时性能好,

另外,开始可使取样频率为1/N,如tt=0:1/2048:6,N=2048.这样显示频率却原频率.

你原程序中第2个频率是150.6*4096/2000=308.4288,取整数+1应是309.

apFFT测相位时,在309点附近任一点测都可以,但有只噪时309点测最好.

你可画出apFFT的振幅谱图和相位谱图,更清楚了.

yangzj

就是说apFFT主要用于测两个同频信号的相位差

但是对于测一个频率未知的正弦信号的起始相位的话我觉得还是有问题:尽管可以通过"已知第N个样点的时间T,从测量值减去相隔T的相位值,即初相位",但这个减去的相位是和正弦信号的频率相关的,而这个精确的频率是并没有得到的.

zhwang554

你对校正很了解,用全相位时移相位差法求出原频率和偏差值,这个精确的频率是可以得到的.只
不过要做2次全相位FFT测相位.

zhwang554

你说"就是说apFFT主要用于测两个同频信号的相位差".

apFFT和FFT测的是样点的即时相位,
某一频率的正弦波的相位是由时间确定的,时间为0时的相位叫初相位,其它时刻的相位就
随时间变化.信号离散化后，每一取样点的时间不同，相位就不同。如一个振幅为1余弦波
，幅值1是另相位，-1是180相位,0值是900或2700,要看后数值关系。apFFT测的是任一样
点时刻的相位,即样点相位。目前流行的相位计测的都是比较相位，被测信号和一个同一
频率参考信号的相位差。

二个同频信号的同一时刻的相位都分别测出来,相位差当然测出了.
可以说apFFT测的是样点的即时相位,可用于测两个同频信号的相位差.

zhwang554

你问"构造信号时，为什么时间序列是这样构造的呢？",


你指是全相位信号的构造?

数字信号处理时有截断和周期展延二个环节,如序列

...x(1) x(2) x(3) X(4) x(5)...

作N=3截断,可以是[x(1) x(2) x(3)],其周期展延是下第1行,可以是[x(2) x(3) x(4)],
其周期展延是下第2行,可以是[x(3) x(4) x(5)],其周期展延是下第3行(原位置展延),

第1行     ...x(1) x(2) x(3) X(1) x(2) x(3)...
第2行     ...x(4) x(2) x(3) X(4) x(2) x(3)...
第3行     ...x(4) x(5) x(3) X(4) x(5) x(3)...

            3X(3)  2X(4)+x(1)  2x(2)+x(5)

这3行相加,形成一个[3X(3)  2X(4)+x(1)  2x(2)+x(5)]周期信号,其1个周期即全相位信
号.
以上是N=3,任何N阶类同,之所以称全相位指对包含x(3)的N=3截断,全部起始相位都考虑
到了.由於所有截断都考虑了,截断引起的Gibbs效应,泄漏可以减轻.
以上是无窗,还可加窗.

我们是为了克服有限数字滤波器的边界效应引入全相位的,构造全相位时间序列
用於FFT是延伸的.

yifeng2030

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)


为什么我加一个n次循环运行后yy2的长度是4096+n，而不是4096呢？

zhwang554

没有看懂你的问题,
下面的程序也可测初相位,避免了循环语句,更简单
你只要对照全相位FFT的方框图,可理解其信号构成过程


NFFT = 4096;
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/NFFT/max(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)

60度

yifeng2030

我的意思是把第一个程序执行一次后，yy2数组的长度是4096+1，而改为：

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



yy2是一个1*4196的数组，而不是1*4096；
还有如果按照全相位的定义，
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 = vecter/NFFT/max(vecter);的意思？
多谢了