[已答复] FFT/APFFT-基于全相位频谱分析的相位差频谱校正法,实际使用时频率如何提取?

[复制链接]
Ljgupup 发表于 2021-12-6 19:42:21
本帖最后由 Ljgupup 于 2021-12-8 08:59 编辑

研究王兆华老师的基于全相位分析的相位差频谱校正法时发现,当采样率为1024时,校正后的频率幅值相位这些参数特别准,但是当我改变采样率时,这个频率怎么校正呢?我试过用谱峰对应的频率加上频偏但是校正后的结果偏差有点大

我的测试结果如下:
采样率fs=1024时校正结果:
f=50.299999,a=0.999988,p=30.000000

采样率fs=512时,相位和幅度都正确,谱峰位置为n=102,频偏102点对应0.6
那么我的频率校正应该怎么计算呢?
fs=512时的频率分辨率为0.5,102*0.5+0.6=51.6,实际为50.3这样计算肯定不对。




我的测试代码如下(参考王兆华老师的代码):

%全相位综合相位差校正法

clear
clc
close all hidden

%准备数据
fs=1024;                %采样频率
nfft=1024;              %FFT计算点数
n=2*nfft-1;             %数据采样点数

%采样间隔定义
t=(-nfft+1):(nfft-1);
%t=0:(2*nfft-2);

%准备数据等间隔采样
x=cos(50.3*2*pi*t/fs+pi/6)+cos(85.3*2*pi*t/fs+pi/3)+cos(121.4*2*pi*t/fs+pi/2);

% =============加窗预处理================
win=hanning(nfft)';

win1=conv(win,win);
win1=win1/sum(win1);
win2=win/sum(win);


%加窗
y1=x.*win1;
y2=x(nfft:2*nfft-1).*win2;

%移位相加
y1a=y1(nfft:end) + [0 y1(1:nfft-1)];


%计算FFT
Y1=fft(y1a, nfft);
Y2=fft(y2, nfft);

%幅度
A1=abs(Y1);
A2=abs(Y2);

%相位
P1=mod(angle(Y1)*180/pi,360);
P2=mod(angle(Y2)*180/pi,360);

%频偏
fe=mod((P2-P1)/180/(1-1/nfft),1);

ae=(A2.^2)./A1*2;

%校正后的值
r=round(50.3);
fre=floor(r)+fe(r+1);
a=ae(r+1);
p=P1(r+1);

sprintf("f=%f,a=%f,p=%f\n", fre, a, p)

22 条回复


zhwang554 发表于 2021-12-8 08:18:56

校正前程序正确,校正部分不正确

本帖最后由 zhwang554 于 2021-12-10 03:45 编辑
Ljgupup 发表于 2021-12-7 16:54
但是问题来了这个频率校正怎么算呢?谱峰位置为n,频率分辨率为f, 是rount(n*f)+频偏吗 ...

修改后测试结果如下:

采样率fs=1024 时校正结果:

f=50.299999,    a=0.999988,    p=30.000000

采样率fs=512 时校正结果:
f=50.300978,    a=1.000002,    p=30.000000

采样率fs=2048 时校正结果:
f=50.299970,   a=1.000057,     p=30.000000


修改后测试代码如下:

clear
clc
close all hidden

%准备数据
fs=1024;                %采样频率
nfft=512;              %FFT计算点数
n=2*nfft-1;             %数据采样点数

%采样间隔定义
t=(-nfft+1):(nfft-1);
%t=0:(2*nfft-2);

%准备数据等间隔采样
x=cos(50.3*2*pi*t/fs+pi/6)+cos(85.3*2*pi*t/fs+pi/3)+cos(121.4*2*pi*t/fs+pi/2);

% =============加窗预处理================
win=hanning(nfft)';

win1=conv(win,win);
win1=win1/sum(win1);
win2=win/sum(win);


%加窗
y1=x.*win1;
y2=x(nfft:2*nfft-1).*win2;

%移位相加
y1a=y1(nfft:end) + [0 y1(1:nfft-1)];


%计算FFT
Y1=fft(y1a, nfft);
Y2=fft(y2, nfft);

%幅度
A1=abs(Y1);
A2=abs(Y2);

%相位
P1=mod(angle(Y1)*180/pi,360);
P2=mod(angle(Y2)*180/pi,360);

%频偏
fe=mod((P2-P1)/180/(1-1/nfft),1);

ae=(A2.^2)./A1*2;

%校正后的值
r=round(50.3*nfft/fs);
fre=(floor(50.3*nfft/fs)+fe(r+1))*fs/nfft;
a=ae(r+1);
p=P1(r+1);

sprintf("f=%f,a=%f,p=%f\n", fre, a, p)



Ljgupup 发表于 2021-12-7 09:30:25
本帖最后由 Ljgupup 于 2021-12-8 09:17 编辑

别沉了!附上图片
apfft1.png

Ljgupup 发表于 2021-12-7 16:44:50
我知道出现这个问题的原因了,问题出在r=round(50.3);这个语句,当采样率为1024时可以这样做,谱峰位置在50这个点,但是当采样率改变之后,其对应的频率分辨率也发生了改变,因此谱峰位置就不是50这个点了。需要重新计算r

Ljgupup 发表于 2021-12-7 16:54:49
Ljgupup 发表于 2021-12-7 16:44
我知道出现这个问题的原因了,问题出在r=round(50.3);这个语句,当采样率为1024时可以这样做,谱峰位置在50 ...

但是问题来了这个频率校正怎么算呢?谱峰位置为n,频率分辨率为f, 是rount(n*f)+频偏吗

Ljgupup 发表于 2021-12-10 15:47:49
zhwang554 发表于 2021-12-8 08:18
修改后测试结果如下:

采样率fs=1024 时校正结果:

谢谢老师,我明白了,今天调试时差法校正也根据你提供的资料调成功了,已经移植到了硬件上去了,准备用实际数据测试一下。感谢您抽空给我解答,我做的产品需要通过测量加速度来精确测量频率,因此您发明的这个算法特别适合。

zhwang554 发表于 2021-12-25 05:26:37

FFT/apFFT校正程序的相位谱和对数振幅谱

本帖最后由 zhwang554 于 2021-12-31 08:01 编辑

aaaa.jpeg

输入信号
exp(j*(2*pi*8.5/N+100*pi/180))+exp(j*(2*pi*18.6/N+100*pi/180))+exp(j*(2*pi*28.7/N+100*pi/180))+exp(j*(2*pi*38.8/N+100*pi/180))
+exp(j*(2*pi*48.9/N+100*pi/180))+exp(j*(2*pi*59.0/N+100*pi/180))+exp(j*(2*pi*69.1/N+100*pi/180))
+exp(j*(2*pi*79.2/N+100*pi/180))+exp(j*(2*pi*89.3N+100*pi/180))+exp(j*(2*pi*99.4/N+100*pi/180))
十个指数信号, 第1个频率8.5,第2个频率18.6,第3个频率28.7,第4个频率38.8, 第5个频率48.9,第6个频率59.0,第7个频率69.1,第8个频率79.2,第9个频率89.3, 第10个频率99.4.   每隔10赫频偏加0.1.     相位 都是100.

图a fft相位谱 p2
X=10,Y=10.3562.        X=20, Y=28.2876,        X=30, Y=46.2164,         X=40,Y=64.1445,          X=50,   Y=82.0723,     
X=60,Y=100,               X=70, Y=117.928,        X=80, Y=135.856,         X=90,Y=153.784,          X=100,  Y=171.716

相邻相位值差值                          p1-p2                                                     频偏
28.2876-10.3562=17.9314           100-10.3562=89..6438=17.9287x5                0.5
46,2164-28.2826=17.9288           100-28.2876=71.7124=17.9031x4                 0.4
64.1445-46.2164=17.9271           100-46.2164=53.7836=17.9279x3                 0.3
82.0723-64.1445=17.9278           100-64.1445=35.8555=17,92775x2               0.2
100-82,0723       =17.9277           100-82.0723=17.9277                                    0.1                                                                                                                                                                   100-100=0                                                      0
117.928-100      =17.928            100-117.928=-17.928                                -0.1
135.856-117.928=17.928             100-135.856=-35.856=-17.928x2                  -0.2
153.784-135.856=17.928             100-153.784=-53.784=17.928x3                   -0.3
171.716-153.784=17.922             100-171.716=-71.716=17.906x4                   -0.4
相邻相位值差值相同                      p1-p2频偏成正比。

图b apfft相位谱 p1
X=10,Y=100       .        X=20, Y=100,              X=30, Y=100,                X=40,Y=100,                 X=50,   Y=100,     
X=60,Y=100,               X=70, Y=100,              X=80, Y=100,                X=90,Y=100,                X=100,  Y=100
表明不信号频偏值多大,apfft相位谱读数不変。

apfft相位用高精度长数据显示如下:
X=   10                             20                            30                            40                            50      
Y=  100.0000000000000   100.0000000000000   100.0000000000001   100.0000000000001  100.0000000000001                          
      60                            70                            80                             90                            100
      100.0000000000000   99.9999999999999    99.9999999999999     99.9999999999999    100.0000000000000
说明apfft相位谱显示精度很高.


图c apfft 对数振幅
     X=10,y=-0.715707     X=20,y=-0.45644     X=30,y=-0.256234     X=40,y=-0.114328     X=50,y=-0.0296292

     X=60,y=-0.0014673    X=70,y=-0.0296292     X=80,y=-0.11429     X=90,y=-0.256082     X=100,y=-.0.455152


图d apfft 对数振幅
      X=10,y=-1.43496          X=20, y=-0.912821      X=30,y=-0.51109       X=40,y=-0.226411     X=50, y=-0.0564929

     X=60,y=1.57374e-07     X=70,y=-0.0564929    X=80,y=-0.226411     X=90,y=-0.51109       X=100, y=-0.912821

图c fft 对数振幅X 2 = 图d apfft 对数振幅

下图是图c和图d的放大对比图。从图中可见图d中apfft的对数振幅谱可达-64.4123db,而fft只能达-30db。这是因为有10个频率成份。旁频泄漏引起的干扰fft比apff严重。
apfftfftooo.jpeg
























aaaa.m

959 Bytes, 下载次数: 9

baa0.m

1.53 KB, 下载次数: 9


Trueternity 发表于 2022-1-12 13:15:27
zhwang554 发表于 2021-12-8 08:18
修改后测试结果如下:

采样率fs=1024 时校正结果:

王老师 能给一个采样点不是从-N开始的全相位时移相位差校正程序借鉴一下  我参照上面程序改了以下 改不太出来 谢谢老师!

Trueternity 发表于 2022-1-13 09:50:16
Ljgupup 发表于 2021-12-7 09:30
别沉了!附上图片

你好 学长 请问一下如果将t改为0:(2*Nfft-2)的话 需要修改校正公式吗 是求频偏的公式会发生变化吗?具体该怎么改呢? 谢谢!!!

zhwang554 发表于 2022-1-14 10:43:40

采样从0开始的全相位时移相位差校正程序

本帖最后由 zhwang554 于 2022-1-23 22:28 编辑
Trueternity 发表于 2022-1-12 13:15
王老师 能给一个采样点不是从-N开始的全相位时移相位差校正程序借鉴一下  我参照上面程序改了以下 改不太 ...

修改后测试结果如下:
f=50.300000,a=0.999985,p=30.000025

clear
clc
close all hidden

%准备数据
fs=512;                %采样频率
nfft=512;              %FFT计算点数
n=2*nfft-1;             %数据采样点数

%采样间隔定义
%t=(-nfft+1):(nfft-1);
t=0:(2*nfft-2);

%准备数据等间隔采样
x=cos(50.3*2*pi*t/fs+pi/6)+cos(85.3*2*pi*t/fs+pi/3)+cos(121.4*2*pi*t/fs+pi/2);

% =============加窗预处理================
win=hanning(nfft)';

win1=conv(win,win);
win1=win1/sum(win1);
win2=win/sum(win);


%加窗
y1=x.*win1;
y2=x(nfft:2*nfft-1).*win2;

%移位相加
y1a=y1(nfft:end) + [0 y1(1:nfft-1)];


%计算FFT
Y1=fft(y1a, nfft);
Y2=fft(y2, nfft);

%幅度
A1=abs(Y1);
A2=abs(Y2);

%相位
P1=mod(angle(Y1)*180/pi,360);
P2=mod(angle(Y2)*180/pi,360);

%频偏
fe=mod((P2-P1)/180/(1-1/nfft),1);

ae=(A2.^2)./A1*2;

%校正后的值
r=round(50.3*nfft/fs);

fre=(floor(50.3*nfft/fs)+fe(r+1))*fs/nfft;

a=ae(r+1);

%P1(r+1)是序列第256样点相位值,利用测到的fre和给定的nfft和fs,反向算出第1样点的相位,即初相位。
%p=mod(P1(r+1)-fre*(nfft-1)/fs*360,360);
%or
p=acosd(cos(P1(r+1)*pi/180 - 2*pi*fre*(nfft-1)/fs));
%or
%p=acosd(cos(angle(Y1(r+1))-2*pi*fre*(nfft-1)/fs))

sprintf("f=%f,a=%f,p=%f\n", fre, a, p)

























补充内容 (2022-3-9 09:59):
用 parzenwin窗 or  bohmanwin窗  代替 hanning 窗 可以改善  频率 振幅相位精度.

zhwang554 发表于 2022-1-15 19:30:42

采样从0开始和从N开始的全相位时移相位差校正程序fft/apfft比较

本帖最后由 zhwang554 于 2022-1-16 07:24 编辑
Ljgupup 发表于 2021-12-10 15:47
谢谢老师,我明白了,今天调试时差法校正也根据你提供的资料调成功了,已经移植到了硬件上去了,准备用实 ...

从0开始测试结果如下:
p   =30.000024674198357          精度小数点后4位
fre =50.299999931327278                               7位
a   = 0.999984534160617                                4位

从-N开始测试结果如下:
p   =30.000000001437005                               8位
fre =  50.29999918532139                               6位
a   =  0.999986326212283                               4位

可见振幅 a 精度都4位相同 频率 fre 精度位6位和7位,基本相同。而相位 p 精度,从0开始4位, 从-N开始8位,-N开始明显改善。这是因为从-N开始的相位值是apfft相位谱直接测量,度高,而从0开始相位值还须校正到第一个样点,校正公中用了频率fre,fre由fft和apfft的相位差测得。fft 的相位精度只有4位,所以总的相位精度就4位差下耒了。
当有多个频率成份时, apfft 的频谱比fft干净. 就是只有一个频率成份如一个余弦信号, 在频谱分析中由於负的鏡像分量的影响,也没有apfft干净.
所以 apfft 频谱的泄漏比 fft 少。所以一篇电子科技大学学位论文评论:
秦磊.滾动轴承振动加速度信号的便携式采集与分析研究[D].电子科技大学, 硕士,2017-4-01,

事实上频谱泄漏引发的谱间干扰问题是影响频谱校正精度的一个重要因素,传统FFT所固有的频谱泄漏效应会很大程序地影响校正法的精度。相比于传统FFT而言全相位FFT具有更优良的抑制频谱泄漏性能,因此可以肯定的是基于全相位FFT谱分析的频滑校正方法比传统频谱校正方法具有更高的校正精度。


fft/apfft实用时,须釆集2N-1个数据存儲后再作频谱分析,所以从-N开始采集能够实现。,


这个fft/apfft方法将fft 和apfft 的特性都用上了,将二者之间的关系也表现出来,特別是振幅校正,是fft振幅谱apfft振幅谱相除。fft/apfft的振幅校正和窗函数无关,别的振幅校正公式都和窗函数有关,而fft/apfft校正方法中的一除将窗函数的影响除掉了, 挺有特色


二种时差校正方法,fft/apfft 和apfft/apfft, 都需要二次FFT変換,许多人都认为硬件实现二次FFT変換计算量大。在许多选用频谱校正的论文中,常见说时差法要二次FFT,一下子就否定了。事实上这二个FFT. 都是实変换,可以合成为一个复数作FFT,FFT是一个复変换,変换后再分开这个方法很有实用意义,和比值法一样,时差法只须一次FFT。

楼主说已经将fft/apfft时差法用於硬件上了, 通过测量加速度来精确测量频率,不知是用二次FFT实现,是用一次FFT实现。


Trueternity 发表于 2022-1-17 14:35:48
zhwang554 发表于 2022-1-14 10:43
采样从0开始的全相位时移相位差校正程序

修改后测试结果如下:

谢谢王老师 我弄明白了 如果采样点为2N-1 则apfft计算出的是第N点的相位 但是我还有一个问题 就是floor(50.3)或是floor(f0)的问题。 如果存在第k次谐波 且同时存在k+0.1次谐波 则使用普通的fft加窗插值算法求解时,k+0.1次谐波会被吞并掉 也就是我无法确定存在k+0.1次谐波所对应的频率;这时候我该如何求解出floor(k+0.1);
换句话说当我面对一个完全未知的信号 在求出频偏后如何确定floor(f0)中的f0具体值是多少呢?
希望老师能够为我解惑!

Trueternity 发表于 2022-1-17 17:14:42
zhwang554 发表于 2022-1-14 10:43
采样从0开始的全相位时移相位差校正程序

修改后测试结果如下:

全相位FFT时移相位差校正方法
https://www.ilovematlab.cn/thread-612948-1-1.html
(出处: MATLAB中文论坛)

老师 我把我有关floor(f0)的问题描述在这个贴子里的品论区里了
在面对未知频率信号的情况下 这个f0是不是通过apfft的幅频特性图中的谱线平均来确定的?
希望老师解答 谢谢老师!

zhwang554 发表于 2022-1-17 21:29:19
本帖最后由 zhwang554 于 2022-1-18 03:25 编辑
Trueternity 发表于 2022-1-17 14:35
谢谢王老师 我弄明白了 如果采样点为2N-1 则apfft计算出的是第N点的相位 但是我还有一个问题 就是floor( ...

程序中的r=round(50.3*nfft/fs) 和floor(50.3*nfft/fs)是振幅特性的峰值处的频率值,不是峰值附近谱线的平均。
如果在一峰值中有二个频率成份,则这个频率分析失效,必须再细分或者加长序列长度以提高谱分辨率。
可以利用apfft相位谱的水平特性耒判别这个峰值是否只有一个频率成份。这个判断方法在一般有噪声的实际信号中也有用。就是在峰值频率对应的相位谱有一个一个像素宽的小平台,若有这个小平台这个峰值是一个频率成分。若是斜线,则是多个频率成分。
参见下贴子:"
http://home.vibunion.com/blog-62061-17761.html全相位水平相位作为判定一个频率成分和信噪比的准则

程序中round和floor 語句在实用程序中是不应该出现的,当时是为了读者易看懂程序没有用搜索峰值的程序,因为在搜索峰值的程序中,有许多判断语句,不易看懂。用判断语句的程序可参见

http://home.vibunion.com/blog-62061-18060.html自动搜索峰值的 apfft/apfft ,fft/apfft 和一位时移4cfft/4cfft校正程序​

Trueternity 发表于 2022-1-18 09:49:34
zhwang554 发表于 2022-1-17 21:29
程序中的r=round(50.3*nfft/fs) 和floor(50.3*nfft/fs)是振幅特性的峰值处的频率值,不是峰值附近谱线的平 ...

我明白了 老师。 我会在设计一个基于labview和matlab混合编程的检测系统中使用到您的算法 用于检测间谐波 真的十分感谢您为我答疑解惑!!! 我也会将您的算法推荐给用于信号处理相关研究方向的同学。
再次感谢!!!

zhwang554 发表于 2022-3-4 05:09:29

不同窗函数的fft/apfft 校正方法校正精度

本帖最后由 zhwang554 于 2022-3-6 21:21 编辑
zhwang554 发表于 2022-1-15 19:30
从0开始测试结果如下:p   =30.000024674198357          精度小数点后4位fre =50.299999931327278        ...


fft/apfft 校正方法利用了fft和apfft 的掁幅谱和相位谱耒进行频率, 相位和振幅的校正,由于相位校正只用apfft 的相位谱读数,不须校
正,所以糈度一直很高, 精度有效位可达12位。但频率校正中用了fft的相位泄㴜誤差大,精度小数点后6-7位。其中振幅校正是fft 振幅谱平方除apfft振幅谱得到,和窗函数无关, 振幅精度为小数点后4位

但matlab测试fft/apfft 校正方法各种窗校正精度值发现.  fft/apfft校正方法校正精度和窗函数有关,测试结果如下
其中parzenwin 和bohmanwin 相位精度达12位, 频率精度达10位,振幅
精度达9位, 精度可以和apfft/apfft-校正法比,fft/apfft须采样2N-1个,apfft/apfft须采样3N-1。这个测试还是初步的,当fs,N,ph1変化精度変化大,也不知是什么原因,但是在有fft频谱校正中振幅和频率精度䏻达到10位,值得注意。

fs = 125; % 采样频率
N = 1024*2;
t = -N+1:N-1;                           % 采样间隔定义

f1 = 50.1234567890123; a1 = 2.1234567890123; ph1 = 100.0000000000351

hanning
fff = 50.123456789015073
aaa = 2.123456790151287
ppp = 1.000000000000339e+02
turkey
fff = 50.12345678907500
aaa = 2.123456803624459
ppp = 1.000000000000338e+02

rectwin
fff =50.123448845879139
aaa = 2.121673434699192
ppp =  1.000000067287741e+02

parzenwin
fff =50.123456789011250
aaa =2.123456789292145
ppp =1.000000000000339e+02

nettaiiwin
fff =50.123456781580806
aaa = 2.123455118540971
ppp =  1.000000000000398e+02

kaiser
fff =50.123449177822948
aaa = 2.121747854113333
ppp =    1.000000061785014e+02

hann
fff =50.123456789073948
aaa =2.123456803389007
ppp =1.000000000000339e+02

hamming
fff =50.123455662637916
aaa =  2.123203575519647
ppp =  1.000000001355033e+02

gausswin
fff =50.123456115432113
aaa =2.123305352725609
ppp = 1.000000000484833e+02

flattopwin
fff =50.123456803267381
aaa = 2.123459994825194
ppp =  1.000000000000560e+02

chevwin
fff = 50.123456816188053
aaa =2.123462900090502
ppp =  1.000000000001128e+02

bohmanwin
fff =50.123456789012614
aaa = 2.123456789083899
ppp =   1.000000000000339e+02

blackmanharris
fff=  50.123456787771673
aaa = 2.123456510054310
ppp =  1.000000000000342e+02

blackman
fff = 50.123483288591757
aaa = 2.123456795130677
ppp =   1.000000000000340e+02

bartlettt
fff = 50.123456800234287
aaa =  2.123459312308074
ppp =   1.000000000000474e+02

barthannwin
fff =50.123456791736182
aaa = 2.123457401484098
ppp = 1.000000000000347e+02

hann(N)=sin(pi*(0:N-1)/(N-1)).^2     
fff =50.123456789055645
aaa =2.123456798758689
ppp = 1.000000000000339e+02

这个测试结果和2010-10-6的测试结果一样,
为什么apfft选用hann窗  http://home.vibunion.com/blog-62061-18152.html
3@blackman                     5.684341886080802e-014    0.09*pi_rad/sample      -100db
5@bohmanwin                 5.684341886080802e-014    0.09*pi_rad/sample       -100db
10@hann                          8.526512829121202e-014   0.06*pi_rad/sample        -100db
13@Parzenwin                5.684341886080802e-014    0.12*pi_rad/sample        -120db

当时测试不同窗的apfft 测相精度。发现现blackman, bohmanwin, hann, Parzenwin 四种窗测相精度达14
当时没有测试测振幅精度,特别是在fft/apfft校正方法中,因为fft/apfft 中测幅是fft振幅谱平方除apfft振蝠谱而得,和窗函数无关。以为两个振愊谱相除将窗函数除掉了。实测显示fft/apfft校正法振幅精和窗函数有关

zhwang554 发表于 2022-3-5 23:41:38

Bohmanwin窗和Parzenwin窗用於楼主贴的fft/apfft程序中

本帖最后由 zhwang554 于 2022-3-15 06:33 编辑
zhwang554 发表于 2021-12-8 08:18
修改后测试结果如下:

采样率fs=1024 时校正结果:

在上贴"不同窗函数的fft/apfft 校正法校正精度"中,parzenwin 和bohmanwin 两个窗精度最好。

将这两种窗用於楼主贴的fft/apfft程序中,测试结果如下(红色位比hanning窗提高位数):
fs=512; nfft=1024*4;      

parzenwin
fm =50.299999999999990
ferror =-7.105427357601002e-15
PH =30.000000000000345
perror =3.446132268436486e-13
AA = 1.499999999996578
Aerror = -3.421707361894732e-12

bohmanwin
fm =50.299999999999919
ferror =-7.815970093361102e-14
PH =30.000000000000338
perror =3.375077994860476e-13
AA =1.499999999997285
Aerror =-2.715383473628208e-12


hanningfm =50.300000000008787
ferror =8.789413641352439e-12
PH =30.000000000000284
perror =2.842170943040401e-13
AA =1.500000000126495
Aerror = 1.264950366675066e-10


Parzenwin窗:  频率精度比hanning窗提高3位,振幅2位,相位0位
Bohmanwin窗:频率精度比hanning窗提高2位,振幅2位,相位0位

看耒,fft/apfft 校正方法的频率和振幅精度的确和窗函数有关。

下面是测试程序,以楼主贴的程序为主,但改为自动搜索峰值(只有一个频率成份)。

clc;clear;clf;close all;
N=1024*4;Fs=512;
t=-N+1:N-1;
f=50.3;
A= 1;
Ph=30;
s=A*cos(2*pi*t*f/Fs+Ph*pi/180);
win=hanning(N)';
win=parzenwin(N)';
win=bohmanwin(N)';
win1=win/sum(win);
win2=conv(win,win);
win2=win2/sum(win2);
s1=s(N:2*N-1);
y1=s1.*win1;
F1=fft(y1,N);
s2=s(1:2*N-1);
y2=s2.*win2;
y2a=y2(N:end)+[0,y2(1:N-1)];
F2=fft(y2a,N);
a1=abs(F1);
[As,k]=max(a1(1:N/2));
P1=mod(angle(F1),2*pi);
P2=mod(angle(F2),2*pi);
dp=P1(k)-P2(k);
dph=mod(dp,2*pi);
if dph>pi
dph=dph-2*pi;
elseif dph<-pi
dph=dph+2*pi;
end   
dk=dph/(pi*(N-1)/N);
AA=abs(F1).^2/abs(F2)*2;
fff=(k+dk-1);
PH=angle(F2);
disp('校正的频率')
fm=[fff]*Fs/N
ferror=fm-f
disp('校正的相位')
PH=mod(PH(k)*180/pi,360)
perror=PH-Ph
disp('校正的振幅')
AA
Aerror=AA-A










zhwang554 发表于 2022-3-12 09:05:54

fft/4cfft频谱校正法

本帖最后由 zhwang554 于 2022-3-30 21:27 编辑
zhwang554 发表于 2022-3-5 23:41
在上贴"不同窗函数的fft/apfft 校正法校正精度"中,parzenwin 和bohmanwin 两个窗精度最好。

将这两种 ...

Apfft 是余弦信号加二个N阶循环卷积窗的fft,更高一阶,4cfft是余弦信号加4个N阶循环巻积窗的fft。
Apfft 和fft结合形成fft/ apfft频谱校正方法。
4cfft和fft结合形成fft/4cfft 频谱校正方法,它仍具有相位不変性,振幅校是fft振幅谱4次和apfft振幅谱相除。频率,柤位,振幅精度很高,三者可同时达到12位以上,它的计算量仍是一次N阶FFT,但輸入数据须4N-3z釆样.
参见 : 4阶移位循环卷积窗ffthttp://home.vibunion.com/blog-62061-18835.html


fft/4cfft频谱校正法测试结果如下:fm =50.299999999999997
ferror =0
PH =0.999999999999427
perrore = -5.726530361016557e-13
AA =1.000000000000000
Aerror =4.440892098500626e-16

clc;clear;clf;close all;
Fs=256;
N=2048*1;
t=0:N*4-2;
t=-2*N+2:2*N-2;
f=50.3;
Ph=1;
A=1;
s=A*cos(2*pi*t*f/Fs+Ph*pi/180);
s=s(1:4*N-3);
win=hanning(N)';
win=sin(pi*(1/2:N-1/2)/N).^2 ;
win=parzenwin(N)';
%win=bohmanwin(N)';
win1=win/sum(win);
win2=conv(win,win);
win4=conv(win2,win2);
win5=win4/sum(win4);
s2=s(1:4*N-3);
s22=s2.*win5;
sb=[0 0 s22(1:N-2)]+s22(N-1:2*N-2)+s22(2*N-1:3*N-2)+[s22(3*N-1:4*N-3) 0];  
F1=fft(sb,N);
a1=abs(F1);
P1=mod(angle(F1),2*pi);
t1=0:N-1;
ss=A*cos(2*pi*t1*f/Fs+Ph*pi/180);  
s1=ss(1:N);
s11=s1.*win1;
F2=fft(s11,N);
a2=abs(F2);
P2=mod(angle(F2),2*pi);
L=fix(N/2)+1;
[A1,k]=max(a1(1:L));
dph=mod(P2(k)-P1(k),2*pi);
if dph>pi
dph=dph-2*pi;
elseif dph<-pi
dph=dph+2*pi;
end   
dk=dph/(pi*(N-1)/N);
fff=k-1+dk;
%AA=(a2(k).^4/a1(k)*8)^(1/3)
AA=(a2(k).^4/a1(k))^(1/3)*2;
PH=P1(k);
fm=[fff]*Fs/N
ferror=fm-f
PH=mod(PH*180/pi,360)
perrore=PH-Ph
AA
Aerror=AA-A



cxzlk 发表于 2022-4-4 16:39:39
Trueternity 发表于 2022-1-18 09:49
我明白了 老师。 我会在设计一个基于labview和matlab混合编程的检测系统中使用到您的算法 用于检测间谐波 ...

同学你好,我想问一下你是怎么处理未知信号频率这个问题的?我最近也在用王老师的全相位这个方法研究间谐波频率,但是不清楚面对一个未知信号怎么去求得他的频率。请解答,谢谢

浅颜轩 发表于 2022-5-12 16:09:03
zhwang554 发表于 2022-1-14 10:43
修改后测试结果如下:
f=50.300000,a=0.999985,p=30.000025

王老师您好,我最近在学习您的全相位fft算法,但是对于算法的校正部分不是很理解,比如这段算法中的y2,Y2,A2,P2这部分,为什么要构造这样的序列?是因为是三角函数而不是指数函数吗?

zhwang554 发表于 2022-5-15 23:06:41

fft和apfft测的相位必须是同一样点的相位,振幅谱必须同样长度

本帖最后由 zhwang554 于 2022-5-18 18:50 编辑
浅颜轩 发表于 2022-5-12 16:09
王老师您好,我最近在学习您的全相位fft算法,但是对于算法的校正部分不是很理解,比如这段算法中的y2,Y2 ...

这个程序是t=0:2*nfft-2的fft/apfft校正法程序,其频率校正公式是程序中语句 fe=mod(P2-P1)/180/(1-1/nfft),1)
其中P2是fft测量的相位,P1是apfft测量的相位,这二个相位必须是一样点的相位

当采样t=0:2*nfft-2时,apfft测得的是中间点第nfft点的相位,所以fft也要测得这个点的相位,由于fft相位谱测得的是序列首位的相位,所以fft序列的nfft采样x(nfft:2*nfft-1):
即程序中语句 :    y2=x(nfft:2*nfft-1).*win2
这样fft测得的P2和apfft测得的P1是同一样点的相位,符合频率校正公式要求。


另外,振幅校正公式是程序中的语句
as=(A2.*2)./A1*2;
即fft振幅谱A2的平方除以apfft振幅谱A1乘2,其中A2和A1必须同样长度才能相除。同时所加窗函数是同一类型,有同样的加权。所以fft序列语句y2=x(nfft:2*nfft-1).*win2
长度为nfft个样点,
而apfft 予处理后是程序中的语句
y1a=y1(nfft:end)+(0 y1(1:nfft-1))
y1a序列长度也是nfft样点长,
as语句中乘2是因为信号是余弦函数而不是指数函数。


所以算法中的y2,Y2,A2,P2这部分,为什么要构造这样的序列?是因为fft apfft 这二个序列测得的相位必须是㒰一样点的,fft和apfft振幅谱必须有相同长度,加同一类型窗函数。







HaoKwok 发表于 2022-5-18 01:49:19
本帖最后由 HaoKwok 于 2022-5-18 12:13 编辑
zhwang554 发表于 2021-12-8 08:18
修改后测试结果如下:

采样率fs=1024 时校正结果:

王老师,我最近也在一直学习您这个apfft
但在我最近使用的实例中我遇到了一个问题就是如下:
我通过给定的x时间序列,然后按照不同的频率w去设计出一个信号y来,然后对y来进行求解w,解出这个信号的频率w是多少。这是我进行程序处理的一个误差图:
为什么我这里的误差这么大呢
附上程序:
clc;
clear;
times = 10;       %验证不同频率w的次数
term = 100;         %测试集频率w的起始点
dterm = 1;     %测试集频率w的间距
real_value_w = (term : dterm : term + dterm*(times-1))';%测试集频率w
get_w = (1:times)';%结果
d_w = (1:times)';%误差
for column= 1 : times
%     column

    %产生原始信号//表达式是固定的
    x_prim=2*pi()*1000/1590:0.0002:2*pi()*1000/1510;        %原始信号时间序列
    simdata = cos(100*x_prim+0*pi()/180);  %原始信号

    %进行插值到2*N-1个点
    N = 2048-1;          %目标插值采样间隔    204

    %插值
    data_x_end = x_prim(end);
    data_x_start = x_prim(1);
    d_u = (data_x_end - data_x_start)/(N - 1);
    x_interp = data_x_start : d_u : data_x_end;
    fs = 1/d_u;
    Y_interp = interp1(x_prim,simdata,x_interp,'spline');
    %插值后的信号以x_interp、Y_interp表示

    %%%%
    %使用apfft解算原始信号的real_value_w的值

    %准备数据
    fs=1000;             %采样频率
    nfft=1024;              %FFT计算点数
    n=2*nfft-1;             %数据采样点数


    f1 = real_value_w(column);

    %         x = cos(2 * pi *x_interp * f1 / fs);
    %cos(real_value_w(column)*x_interp+0*pi()/180);
    x = Y_interp;

    % =============加窗预处理================
    win=hanning(nfft)';

    win1=conv(win,win);
    win1=win1/sum(win1);
    win2=win/sum(win);
    %加窗
    y1=x.*win1;
    y2=x(nfft:2*nfft-1).*win2;

    %移位相加
    y1a=y1(nfft:end) + [0 y1(1:nfft-1)];

    %计算FFT
    Y1=fft(y1a, nfft);
    Y2=fft(y2, nfft);

    %幅度
    A1=abs(Y1);
    A2=abs(Y2);

    %相位
    P1=mod(angle(Y1)*180/pi,360);
    P2=mod(angle(Y2)*180/pi,360);

    %频偏
    fe=mod((P2-P1)/180/(1-1/nfft),1);

    ae=(A2.^2)./A1*2;

    %校正后的值
    r=round(f1*nfft/fs);
    fre=(floor(f1*nfft/fs)+fe(r+1))*fs/nfft;
    a=ae(r+1);
    p=P1(r+1);
    get_w(column)=fre;
    d_w(column)=real_value_w(column)-fre;
end
real_value_w    %未知信号的频率
get_w           %求解出的频率
d_w             %误差




zhwang554 发表于 前天 22:59
本帖最后由 zhwang554 于 2022-5-27 05:29 编辑
HaoKwok 发表于 2022-5-18 01:49
王老师,我最近也在一直学习您这个apfft
但在我最近使用的实例中我遇到了一个问题就是如下:
我通过给定的 ...

如何证明你构造的信号的采样频率fs=1000  ?


对你的信号构成过程看不懂。
频率100时的2048个采样的x=Y-interp的图形是附图上图,在整个2048采样序列中约有三个周期多一点,这个信号你先产生1047个采样的 x-prim 和simdate ,再通过插值产生Y-interp, 然后你认为构成一个2048个采样的fs=1000  的100赫余弦信号,

一个正常的100赫余弦信号耍在2048采样中产生三个周期多一点,取样频率约为60000 多。附图下图是cos(2*pi*100*(0:2047)/61400-38.491*pi/180)的图形,二者很相近。

所以你如何证明你构造的信号的采样频率fs=1000? nnnn.jpg

您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表