|
本帖最后由 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-10 03:45 编辑 Ljgupup 发表于 2021-12-7 16:54 修改后测试结果如下： 采样率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-8 09:17 编辑 别沉了！附上图片 |
| 我知道出现这个问题的原因了，问题出在r=round(50.3);这个语句，当采样率为1024时可以这样做，谱峰位置在50这个点，但是当采样率改变之后，其对应的频率分辨率也发生了改变，因此谱峰位置就不是50这个点了。需要重新计算r |
Ljgupup 发表于 2021-12-7 16:44 但是问题来了这个频率校正怎么算呢？谱峰位置为n，频率分辨率为f, 是rount(n*f)+频偏吗 |
zhwang554 发表于 2021-12-8 08:18 谢谢老师，我明白了，今天调试时差法校正也根据你提供的资料调成功了，已经移植到了硬件上去了，准备用实际数据测试一下。感谢您抽空给我解答，我做的产品需要通过测量加速度来精确测量频率，因此您发明的这个算法特别适合。 |
本帖最后由 zhwang554 于 2021-12-31 08:01 编辑 
输入信号 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严重。 
|
zhwang554 发表于 2021-12-8 08:18 王老师 能给一个采样点不是从-N开始的全相位时移相位差校正程序借鉴一下 我参照上面程序改了以下 改不太出来 谢谢老师！ |
Ljgupup 发表于 2021-12-7 09:30 你好 学长 请问一下如果将t改为0：(2*Nfft-2)的话 需要修改校正公式吗 是求频偏的公式会发生变化吗？具体该怎么改呢？ 谢谢！！！ |
本帖最后由 zhwang554 于 2022-1-23 22:28 编辑 Trueternity 发表于 2022-1-12 13:15 修改后测试结果如下： 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-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实现。 |
zhwang554 发表于 2022-1-14 10:43 谢谢王老师 我弄明白了 如果采样点为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具体值是多少呢？ 希望老师能够为我解惑！ |
zhwang554 发表于 2022-1-14 10:43 全相位FFT时移相位差校正方法 https://www.ilovematlab.cn/thread-612948-1-1.html (出处: MATLAB中文论坛) 老师 我把我有关floor（f0）的问题描述在这个贴子里的品论区里了 在面对未知频率信号的情况下 这个f0是不是通过apfft的幅频特性图中的谱线平均来确定的？ 希望老师解答 谢谢老师！ |
本帖最后由 zhwang554 于 2022-1-18 03:25 编辑 Trueternity 发表于 2022-1-17 14:35 程序中的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校正程序 |
zhwang554 发表于 2022-1-17 21:29 我明白了 老师。 我会在设计一个基于labview和matlab混合编程的检测系统中使用到您的算法 用于检测间谐波 真的十分感谢您为我答疑解惑！！！ 我也会将您的算法推荐给用于信号处理相关研究方向的同学。 再次感谢！！！ |
本帖最后由 zhwang554 于 2022-3-6 21:21 编辑 zhwang554 发表于 2022-1-15 19:30 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-15 06:33 编辑 zhwang554 发表于 2021-12-8 08:18 在上贴＂不同窗函数的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-30 21:27 编辑 zhwang554 发表于 2022-3-5 23:41 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 |
Trueternity 发表于 2022-1-18 09:49 同学你好，我想问一下你是怎么处理未知信号频率这个问题的？我最近也在用王老师的全相位这个方法研究间谐波频率，但是不清楚面对一个未知信号怎么去求得他的频率。请解答，谢谢 |
zhwang554 发表于 2022-1-14 10:43 王老师您好，我最近在学习您的全相位fft算法，但是对于算法的校正部分不是很理解，比如这段算法中的y2,Y2,A2,P2这部分，为什么要构造这样的序列？是因为是三角函数而不是指数函数吗？ |
本帖最后由 zhwang554 于 2022-5-18 18:50 编辑 浅颜轩 发表于 2022-5-12 16:09 这个程序是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 12:13 编辑 zhwang554 发表于 2021-12-8 08:18 王老师，我最近也在一直学习您这个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 于 2022-5-27 05:29 编辑 HaoKwok 发表于 2022-5-18 01:49 如何证明你构造的信号的采样频率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? 
|
微博
QQ好友和群
来自 5#
