FFT加窗插值修正后，求解多次谐波频率，幅值，相位的精度杠杠的（暂时不开源）

eric2013

此贴算法暂时没有开源打算（2022-07-18添加）

以工频信号为例，如果波形是固定的50Hz，我们可以直接设置采样率是这个周期的整数倍，并进行同步采样，求出来的数据精度就比较高，但实际频率是50Hz左右，从会导致FFT求解出来的频率，幅值，相位都有偏差。

所以就有了各种加窗处理（其实就是加权），FFT处理完毕后再做插值修正（目的是解决栅栏效应），准确率高了很多。



设置信号由10次谐波组成，基波49.9Hz。
谐波幅值 =[11,2,1,0.5,0.3,0.1,0.05,0.05,0.05,0.05];
谐波初始相位=[0.9,0.8,0.7,0.6,0.5,0.3,0.2,0.2,0.2,0.2];


信号先做汉宁窗，FFT后再做插值，10次谐波的频率，赋值和相位如下：
   1     49.901010   11.000616  0.899702
   2     99.768131   1.999112    0.808405
   3     149.673843  0.999671   0.706992
   4     199.562452  0.499779  0.610078
   5     249.454945  0.299804  0.512082
   6     299.301802  0.099861  0.327334
   7     349.223025  0.049967  0.221232
   8     399.165802  0.050010  0.209194
   9     449.065492  0.050015  0.209054
  10     498.939165  0.049954  0.215565

艾那的小强

只知道加窗，插值是用来干什么的？是用来提高精度么？

eric2013

艾那的小强 发表于 2020-5-17 23:17
只知道加窗，插值是用来干什么的？是用来提高精度么？

加窗是解决频谱泄露问题，插值是解决栅栏效应。

missfox

继续坐等白嫖

一次嫖一次爽，一直嫖一直爽

艾那的小强

eric2013 发表于 2020-5-18 01:35
加窗是解决频谱泄露问题，插值是解决栅栏效应。

硬汉哥。能请教一下，怎么对FFT结果插值吗？

lvehe

对音频信号进行FFT时存在以下问题，请问如何解决。
音频采样频率为48000Hz，FFT为1024线，48000/1024=46.875Hz。
不是整数，如果是50次就完美了，因为要求1秒出50个频谱。
一种方法是采用升采样的方法，比如实际获得960点数据即开始FFT计算。
但是FFT算法必须是2的N次方，而960不是，所以计算之前要将数据扩展为1024个。
请问如何对数据进行升采样操作？数据扩充后，FFT计算得到的频率、幅值、能量等是否要修正？如何修正？

eric2013

lvehe 发表于 2020-5-20 11:16
对音频信号进行FFT时存在以下问题，请问如何解决。
音频采样频率为48000Hz，FFT为1024线，48000/1024=46.8 ...

不用修正，补零后前后效果。

FFT补零计算的前后效果
http://www.armbbs.cn/forum.php?m ... 7804&fromuid=58
(出处: 硬汉嵌入式论坛)

eric2013

艾那的小强 发表于 2020-5-20 09:56
硬汉哥。能请教一下，怎么对FFT结果插值吗？

下次DSP教程更新。

艾那的小强

eric2013 发表于 2020-5-20 11:34
下次DSP教程更新。

哦。谢谢硬汉哥。

艾那的小强

eric2013 发表于 2020-5-20 11:34
下次DSP教程更新。

请问一下，硬汉哥。你楼主位的例子后面会给出来吗？想学习一下。

eric2013

艾那的小强 发表于 2020-5-21 18:07
请问一下，硬汉哥。你楼主位的例子后面会给出来吗？想学习一下。

会。下期DSP先发布matlab操作，后面再发布板子操作。

艾那的小强

eric2013 发表于 2020-5-22 07:58
会。下期DSP先发布matlab操作，后面再发布板子操作。

期待。谢谢硬汉哥。

xiaosir

硬汉哥,我语音信号FFT->填充滤波->IFFT后有杂音,间断的哒声,周期跟每次处理的长度一致,是不是也可以这样处理?

eric2013

xiaosir 发表于 2020-5-23 11:44
硬汉哥,我语音信号FFT->填充滤波->IFFT后有杂音,间断的哒声,周期跟每次处理的长度一致,是不是也可以这样处 ...

不行，这个是为了更好的求解谐波的幅值和相位。

艾那的小强

硬汉哥，我现在用了你说的，加了汉宁窗，也做了全相位处理，现在算手动生成的数据（比如我在程序里面生成一组正弦波数据），相位和频率均能计算正确。

但我用ADC采集信号发生器发出的正弦波，却总是测不准。是ADC采集这里有什么特殊要求吗？

我用的是你例子里面的定时器触发ADC采集，DMA中断的方式。在DMA完成中断中关掉定时器。然后计算FFT（在主循环中计算FFT)。计算完了再打开定时器。
能给点指点吗？谢谢了。

eric2013

艾那的小强 发表于 2020-5-25 14:26
硬汉哥，我现在用了你说的，加了汉宁窗，也做了全相位处理，现在算手动生成的数据（比如我在程序里面生成一 ...

等我发布了，参考即可。

艾那的小强

eric2013 发表于 2020-5-25 15:44
等我发布了，参考即可。

硬汉哥，我想调皮一下。
能不能开通个付费尝鲜服务。
我被这个问题卡了两周了。实在是着急了。

eric2013

艾那的小强 发表于 2020-5-25 17:31
硬汉哥，我想调皮一下。
能不能开通个付费尝鲜服务。
我被这个问题卡了两周了。实在是着急了。

没有这种服务，只能等。

艾那的小强

eric2013 发表于 2020-5-25 17:32
没有这种服务，只能等。

嗯嗯。好的。那我再查查。期待。。

艾那的小强

顶一下，硬汉哥别忙忘了这事。

haci

版主，请教下，如果采样频率不是（幅值最大频率）的整数倍，而且不是50.05倍这种很小的，比如说50.5倍，通过加窗和插值可以解决吗？
我现在加了ntuuall窗，再用的时域平移相位差校正法，效果还是不理想，具体表现为（采样频率不是整数倍时），计算相位与实际相位差90度，
（采样频率是整数倍时），计算相位与实际相位无差别，我现在的做法是：计算出频率后(频率计算的比较准)，根据当前频率实时调整采样频率，这样才能得到比较好的效果
PS:信号频率是变化的。
请教下有没有更好的办法

eric2013

haci 发表于 2020-7-1 22:18
版主，请教下，如果采样频率不是（幅值最大频率）的整数倍，而且不是50.05倍这种很小的，比如说50.5倍，通 ...

此贴就是解决这个问题而生，待我后面腾出精力了，分享相关教程和Demo到第2版DSP教程里面。

migsoft

xiaosir 发表于 2020-5-23 11:44
硬汉哥,我语音信号FFT->填充滤波->IFFT后有杂音,间断的哒声,周期跟每次处理的长度一致,是不是也可以这样处 ...

FFT和IFFT要进行重叠处理

ou513

佬大这个demo例程能发布一下，参考参考

eric2013

ou513 发表于 2020-8-14 21:28
佬大这个demo例程能发布一下，参考参考

下次DSP教程更新时发布matlab版，后面发布单片机版。

zhm19840622

非常期待硬汉MATLAB版例程！！！！！！！！！！！！！

eric2013

zhm19840622 发表于 2020-8-24 20:59
非常期待硬汉MATLAB版例程！！！！！！！！！！！！！

好，正确早点发布。

zhm19840622

eric2013 发表于 2020-8-25 08:56
好，正确早点发布。

牛！！！！！！！！！

713042507

eric2013 发表于 2020-8-25 08:56
好，正确早点发布。

跪求硬汉哥早点发布，等着呢

713042507

硬汉哥，小白请教一下，
我现在需要做两个正弦波的相位差，
x=SIN(2*PI* 64)，初始相位为0；
算4096个点的FFT,分辨率为1HZ，然后求出幅值，
通过比较法；找到幅值最大的那个点对应的频率。
然后在将频率带入到FFT的输出数组算出来相位为-87度



代码如下：
      for (i=0;i<= FFT_LENGTH ; i++)      
      {
          fft_inputbuf [2*i  ] = data;
          fft_inputbuf [2*i+1] = 0;
      }
      arm_cfft_radix4_f32(&scfft,fft_inputbuf); //FFT计算（基4）

      arm_cmplx_mag_f32(fft_inputbuf,fft_outputbuf_amp,FFT_LENGTH); //算出幅值

      dc_vaule = fft_outputbuf_amp[0] /4096;//算出直流分量
   
      for (i=0;i<= FFT_LENGTH/2 ; i++)
      {
        fft_outputbuf_amp = fft_outputbuf_amp[i+1] ;  //去掉直流分量，去前2047个点
      }

      arm_max_f32(fft_outputbuf_amp,2047,&max_seq_value,&mean_value_seq); //在2047个点中找到幅值最大的那个点
   
      phase = atan2(fft_inputbuf[(mean_value_seq+1)*2+1], fft_inputbuf[(mean_value_seq+1)*2]);  //算出幅值最大的点的相位，

713042507

艾那的小强 发表于 2020-5-25 14:26
硬汉哥，我现在用了你说的，加了汉宁窗，也做了全相位处理，现在算手动生成的数据（比如我在程序里面生成一 ...

小哥，能借鉴一下你的FFT求解相位的代码吗？我现在频率求得很准，但相位算不对

滴滴滴

谢谢！！！！

dufanbao

期待 顶一下

嵌入式小白

哥 MATLAB程序有了吗？参考一下

嵌入式小白

713042507 发表于 2020-9-9 20:54
小哥，能借鉴一下你的FFT求解相位的代码吗？我现在频率求得很准，但相位算不对

我也想要哥

嵌入式小白

为什么我加窗之后 栅栏效应更明显了？ 代码如下：
%% 窗函数测试
clc
close all

Ts=12800;
Fs=1/Ts;
t=1/12800:1/12800:0.02;
yt = sin(2*pi*50*t)+ 0.1*sin(3*2*pi*50*t)+ 0.5*sin(5*2*pi*50*t);%原始电网信号
fft_out1 = abs(fft(yt));
figure
subplot(211)
plot( yt)
title('原始信号')
subplot(212)
stem(fft_out1)
title('原始信号频谱')

win = hanning(length(t));%调用汉宁窗。
yt1 = yt.*win';
fft_out2=abs(fft(2*yt1));
figure
subplot(211)
plot(yt1)
title('加窗信号')
subplot(212)
stem(fft_out2)
title('加窗信号频谱');

eric2013

嵌入式小白 发表于 2020-11-6 12:27
为什么我加窗之后 栅栏效应更明显了？ 代码如下：
%% 窗函数测试
clc

过段时间时间开源发布后，参考即可。

嵌入式小白

eric2013 发表于 2020-11-6 13:24
过段时间时间开源发布后，参考即可。

哥 啥时候出啊 好期待啊

yaoyfuan

期待期待！

yaoyfuan

硬汉哥，啥时候发布matlab版呀，最近一直在学，看不明白