求助：离散信号全相位fft求相位

qiuyun0214

已知一随机离散信号，周期为16，采样频率6k，采了2000个周期共2000*16个数据，现在用fft和apfft分别用五种方法求相位，相位值都有差别，请问哪个最为精确？
程序如下：
clear all;
clc;
Fs=6000;
%“data.mat”内存放的是2000*16矩阵，即以16个点为周期采样频率为6k采样的数据
load data.mat;
%--fft1 直接对平均后的周期信号作16点fft（默认矩形窗）
v=mean(m);
v1=v-mean(v); %去均值
y1=fft(v1,16);      
p1=angle(y1(2));%幅值最大处相位
fprintf('p1=%f\n',p1);

%--fft2 用hanning窗对平均后的周期信号作16点fft
[A P]=FFT_window(v1,16);
rw=round(375/Fs*16)+1;%谱峰值所处的谱线
p2=P(rw);
fprintf('p2=%f\n',p2);

%--fft3 对长为2000*16的原信号进行fft
v3=m';v3=v3(:);v3=v3-mean(v3);v3=v3';
[A PP]=FFT_window(v3,length(v3));
rw=round(375/Fs*length(v3))+1;%谱峰值所处的谱线
p3=PP(rw);
fprintf('p3=%f\n',p3);

%--apfft1 平均后对周期信号作延拓
vv=mean(m);vv(17)=vv(1);vv=vv-mean(vv);%用一个数补上末尾
V1(17:17*2-1)=vv(1:17);
V1(1:16)=vv(1:16);

[A AP1]=apFFT_analysis(V1,17);
rw=round(375/Fs*17)+1;%谱峰值所处的谱线
ap1=AP1(rw);
fprintf('ap1=%f\n',ap1);

%--apfft2 直接apfft
vv2=m';vv2=vv2(:);vv2(length(vv2)+1)=vv2(1);vv2=vv2-mean(vv2);
NN=length(vv2);

V2(NN:NN*2-1)=vv2(1:NN);
V2(1:NN-1)=vv2(1:NN-1);
[A AP2]=apFFT_analysis(V2,NN);
rw=round(375/Fs*NN)+1;%谱峰值所处的谱线
ap2=AP2(rw);
fprintf('ap2=%f\n',ap2);


function  [a,p]= apFFT_analysis(v,NN)
    win =  hanning(NN)';
    winn=conv(win,win);%apFFT需要卷积窗
    winap=winn/sum(winn);%窗归一化
    Y=v.*winap;
    YY=Y(NN:end)+[0 Y(1:NN-1)];%构成长N的FFT输入数据
    YY_fft=fft(YY,NN);
    a=abs(YY_fft);%apFFT的振幅谱      
    p=mod(phase(YY_fft)*180/pi,360);%apFFT的相位谱
    p=p*pi/180;

function  [a,p]= FFT_window(v,NN)
    win =  hanning(NN)';
    win1 = win/sum(win);%窗归一化
    Y= v.*win1;
    Y_fft=fft(Y,NN);   
    a=abs(Y_fft);
    p=mod(phase(Y_fft)*180/pi,360);
    p=p*pi/180;

得出来的相位值如下：
p1=0.292931
p2=0.293172
p3=0.292881
ap1=0.300581
ap2=0.292915
        

qiuyun0214

将“data.dat”改成“data.mat"即可。其中p1和p2的区别是因为fft默认为矩形窗，即p1是用矩形窗求出的，而p2是用hanning窗求出的，这两种方法是将周期信号平均后成为一个周期来求的；第三种是直接对信号进行fft，fft计算数据量大，速度很慢；其中有一个问题是：全相位fft需要对数据进行预处理，即进行周期延拓，而在我的离散信号中第一点与最后一点是不一样的，所以我第一步将数据的第一点补充到最后一点形成2000*16+1个数据，然后再进行周期延拓就能产生奇数个数据，不知道这样延拓是不是正确的？
究竟是平均了再求比较好，还是直接求比较好？因为结果值确实都有差异，有点弄糊涂了！望指教！

qiuyun0214

其中375是基频，即采样频率Fs/周期N=375

zhwang554

回复 2 # qiuyun0214 的帖子

全相位fft只能处理稳定确定信号,即三角函数信号不能处理随机信号
你说:”全相位fft需要对数据进行预处理，即进行周期延拓”, 信号周期延拓后再作全相位处理,得到的还是原信号.所以你的全相位fft和fft一样,不会有区别.
全相位fft数据不能通过周期延拓得到,必连续采样得到
对整周期采样信号,apfft和fft是相同的

qiuyun0214

回复 4 # zhwang554 的帖子

非常感谢老师的指导！
原来全相位fft不能处理随机信号，难怪我在看到所有的例子时都是对标准三角函数信号进行处理的，原来还是自己理解不够深刻。导致将例子套用时，就出现了周期延拓的问题。
    再请教老师一个稍微前沿一点的问题，目前有没有办法能够解决同频噪声的问题？在信噪比较高时，同频干扰基本可以忽略；而在信噪比较低的情况下，同频噪声会产生很严重的破坏性干扰，在软件方法上有没有办法能够减弱这种同频干扰呢？

wwllzz

我在计算用全相位预处理再傅立叶解析, 算相位, 要想真正实现高分辨率还真不容易. 处理器的字长受限制, atan2在某些点附近本身的运算误差很大, 哪怕双精度浮点数.  我是避开某些相位 多算几个平均. 明天看看结果如何

wwllzz

我采用全相位预处理, 采样点和周期的的误差控制在0.5点以内, 我算下来主要是结果跳动难以忍受 整周波采样附近很稳定10万分之一, 差0.5点能达到万分之一就不错了. 考虑下加窗吧 再看看结果如何

qiuyun0214

回复 7 # wwllzz 的帖子

你计算的数据是什么数据？我发上去的还是你自己的？

wwllzz

我自己的数据 加窗, 也许是我窗加的不对 效果更差 稳定是大问题 稳定的误差才能补偿

wwllzz

用相位差计算频率, 最差条件(差0.5采样点)需要很长时间才能稳定到万分之0.5以内(50Hz 12秒),  非整周波采样条件下, 如何令相位快速准确稳定的估算出来.  还是我的算法程序写的有问题

qiuyun0214

回复 10 # wwllzz 的帖子

我不太明白你的意思，你的“非整周波采样”是指什么？你说的最差条件差0.5采样点计算出来的频率，这种误差的大小应该和采样频率也有关。而且地板老师说了，全相位只适合处理稳定确定信号，不能处理随机信号，你用全相位计算和fft计算差别就不大了，那么分辨率肯定是不高的，你可以考虑离散频谱校正，论坛中有个帖子是关于相位差法之改变窗长法来进行离散频谱校正，你可以参考参考。

wwllzz

回复 11 # qiuyun0214 的帖子

现阶段我只考虑相位校正, 我处理的算稳定确定信号吧, 就是电网的50Hz, 这频率的正常允许波动范围是+/-0.5Hz, 说分辨率是在实际条件下, 加上稳定的模拟信号, 在多次测量下的频率或相位差应该是稳定的值, 而不是忽大忽小的跳变量. 比如测量A相电压的相位和B相电压的相位应该稳定的差120度对应的时间差(相位差), 实际做出来的结果则出乎理论的意料之外. 我希望做到这差值在1uS以内, 万分之0.5

qiuyun0214

回复 12 # wwllzz 的帖子

这我就搞不懂了啊，既然你频率有正常波动范围，那么从高精度上来讲，相位也会有偏移，那A相电压相位和B相电压相位的相位差当然就有可能是呈现忽大忽小的跳变趋势了啊？如果确实是相位差变化了，那你用什么方法做出来的结果当然都是跳变的。那你不同情况下测量的信号有噪声么？我不清楚你这个领域信号情况。

zhwang554

   不清楚为什么你的数据整周波采样附近很稳定10万分之一, 差0.5点能达到万分之一就不错了，

   你可能测的是A相电压的相位和B相电压的相位, 应该稳定的差120度,

   下表一和表二为利用本论坛日志

fft/apfft谱分析实例一----测3相交流电的相位2009-04-12 06:21

http://forum.vibunion.com/home-space-uid-62061-do-blog-id-17736.html

中的三相电数据data1.txt, 是一个3X5305采样数据.

    在不同起点n,测量三相频率f1,f2,,f3,及相位差p1-p2 , p2-p3 , p3-p1 如表一表二(N=1024阶apfft/apfft校正法)

           表一   不同起点三相电频率

n                    f1                               f2                               f3

0     50.036090281941    50.0367483793498    50.0195615359034

2     50.0361232479126   50.0369093887453    50.0195188240093

20    50.0363594446624   50.0382729216715    50.0191757617334

100   50.0360598760094   50.0419632458548    50.0184712742239

200   50.0331435117517   50.039995942718     50.0188442345439

500   50.0288750332103   50.0157828346402    50.018093427237

1000  50.0329828640265   50.0281936615661    50.0342013546007

2000  50.029599010895    50.0314038332878    50.0379686574494

           表二   不同起点三相电相位差

n              p1-p2                            p2-p3                       p3-p1    

0     -240.882268366594   122.860738844139    118.021529522455

2     -240.872423581266   122.848525294691    118.023898286576

20    -240.787779993835   122.745428648124    118.042351345711

100   -240.511161141149   122.462156387821    118.049004753329

200   -240.407451719457   122.551489962266    117.85596175719

500   -241.163338935935   123.826834717937    117.336504217997

1000  -240.93985893497    124.090175960394    116.849682974575

2000  -240.627271305199   123.779448769737    116.847822535462

    由表一表二可见,同一点三相频率f1,f2,,f3不相同, 如f1和f2差0.00x, 所以各点p1-p2相位差不会相同.只有f1和f2严格同频, 各点p1-p2相位差才相同

     不同点同一频率f不相同, 如f1,差0.00x

     所以你测的是A相电压的相位和B相电压的相位, 应该稳定的差120度, 我估计A相和B相的频率不严格同频,所以相位差精度能达到万分之一就不错了.

         A相和B相的频率不严格同频可能是测量噪声引起, 但不明白为什么你测的整周波采样附近很稳定10万分之一,

wwllzz

感谢老师能关注我的留言. 我是用实际的实验设备加模拟信号, 再AD采样计算的, 不是纯数学模拟. atan函数计算任意角度实际运算误差就很大, 我原来的算法也有些误区, 就是计算信号的任意初相位值的实际计算误差可能差很多. 比如atan算初相位接近90度附近精度难以提高, 后来采用计算小的相位差用asin函数来求,  得出的结果稳定性好些.  我会多留意下两相信号采样不严格同频率的问题.  实际的系统环境下远比数学模拟复杂, 能尽快抓到问题要害的好, 有新的体会我继续和大家分享