fft测量频率不准问题

sssuuu · 发表于 2025-3-19 17:25:28

本帖最后由 sssuuu 于 2025-3-19 17:27 编辑

使用640ksps，1024个点采样，分辨率为640Hz，为什么采样5kHz时是准的，采样10kHz、15kHz。。。50kHz。。。时都不准，我看了频谱，基波索引都是提前了，就是测得小了。

例如采样10kHz，测得9375Hz
采样50kHz，测得48125Hz，这也差的太多了吧。
是fft运算的问题吗

这里是我fft代码，用的是硬汉哥的，文件上传不上来：

/*
*********************************************************************************************************
* 函数名: Int_FFT_TAB
* 功能说明: 正弦和余弦表
* 形参: FFT点数
* 返回值: 无
*********************************************************************************************************
*/
#if (MAX_FFT_N != 8192) && (MAX_FFT_N != 16384) && (MAX_FFT_N != 32768)
float32_t costab[MAX_FFT_N/2];
float32_t sintab[MAX_FFT_N/2];
void InitTableFFT(uint32_t n)
{
uint32_t i;

/* 正常使用下面获取cos和sin值 */
#if 1
for (i = 0; i < n/2; i ++ )
{
sintab[ i ]=  sin( 2 * PI * i / MAX_FFT_N );
costab[ i ]=  cos( 2 * PI * i / MAX_FFT_N );
}

/* 打印出来是方便整理cos值和sin值数组，将其放到内部Flash，从而节省RAM */
#else
Set_Current_USART(USART1_IDX);
printf("const float32_t sintab[%d] = {\r\n", n);
for (i = 0; i < n/2; i ++ )
{
sintab[ i ]=  sin( 2 * PI * i / MAX_FFT_N );
printf("%.11ff,\r\n", sintab[ i ]);
}
printf("};\r\n");

printf("const float32_t costab[%d] = {\r\n", n);
for (i = 0; i < n/2; i ++ )
{
sintab[ i ]=  cos( 2 * PI * i / MAX_FFT_N );
printf("%.11ff,\r\n", sintab[ i ]);
}
printf("};\r\n");
#endif
}
#endif

/*
*********************************************************************************************************
* 函数名: cfft
* 功能说明: 对输入的复数组进行快速傅里叶变换（FFT）
* 形参: *FFT_Input  指向复数数组的指针，数组中实部和虚部交替存储，即 [real0, imag0, real1, imag1, ...]
*          FFT_N 表示点数，必须是2的幂次方
* 返回值: 无
*********************************************************************************************************
*/
void cfft(float32_t *FFT_Input, uint32_t FFT_N )
{
float32_t TempReal1, TempImag1, TempReal2, TempImag2;
uint32_t k,i,j,z;
uint32_t Butterfly_NoPerColumn;       /* 每级蝶形的蝶形组数 */
uint32_t Butterfly_NoOfGroup; /* 蝶形组的第几组 */
uint32_t Butterfly_NoPerGroup; /* 蝶形组的第几个蝶形 */
uint32_t ButterflyIndex1,ButterflyIndex2,P,J;
uint32_t L;
uint32_t M;

z=FFT_N/2;                   /* 变址运算，即把自然顺序变成倒位序，采用雷德算法 */

for(i=0,j=0;i<FFT_N-1;i++)
{
/*
   如果i<j,即进行变址 i=j说明是它本身，i>j说明前面已经变换过了，不许再变化，注意这里一般是实数有虚数部分设置成结合体
*/
if(i<j)
{
TempReal1  = FFT_Input[j*2];
FFT_Input[j*2]= FFT_Input[i*2];
FFT_Input[i*2]= TempReal1;
}

k=z;                        /*求j的下一个倒位序 */

while(k<=j)                   /* 如果k<=j,表示j的最高位为1 */
{
j=j-k;                     /* 把最高位变成0 */
k=k/2;                     /* k/2，比较次高位，依次类推，逐个比较，直到某个位为0，通过下面那句j=j+k使其变为1 */
}

j=j+k;                   /* 求下一个反序号，如果是0，则把0改为1 */
}

/* 第L级蝶形(M)第Butterfly_NoOfGroup组(Butterfly_NoPerColumn)第J个蝶形(Butterfly_NoPerGroup)****** */
/* 蝶形的组数以2的倍数递减Butterfly_NoPerColumn，每组中蝶形的个数以2的倍数递增Butterfly_NoPerGroup */
/* 在计算蝶形时，每L列的蝶形组数,一共有M列，每组蝶形中蝶形的个数,蝶形的阶数(0,1,2.....M-1) */
Butterfly_NoPerColumn = FFT_N;
Butterfly_NoPerGroup = 1;
M = (uint32_t)log2(FFT_N);
for ( L = 0;L < M; L++ )
{
Butterfly_NoPerColumn >>= 1; /* 蝶形组数假如N=8，则(4,2,1) */

/* 第L级蝶形第Butterfly_NoOfGroup组（0,1，....Butterfly_NoOfGroup-1）*/
for ( Butterfly_NoOfGroup = 0;Butterfly_NoOfGroup < Butterfly_NoPerColumn;Butterfly_NoOfGroup++ )
{
for ( J = 0;J < Butterfly_NoPerGroup;J ++ )       /* 第 Butterfly_NoOfGroup 组中的第J个 */
{          /* 第 ButterflyIndex1 和第 ButterflyIndex2 个元素作蝶形运算,WNC */
ButterflyIndex1 = ( ( Butterfly_NoOfGroup * Butterfly_NoPerGroup ) << 1 ) + J;/* (0,2,4,6)(0,1,4,5)(0,1,2,3) */
ButterflyIndex2 = ButterflyIndex1 + Butterfly_NoPerGroup;/* 两个要做蝶形运算的数相距Butterfly_NoPerGroup (ge=1,2,4) */
P = J * Butterfly_NoPerColumn; /* 这里相当于P=J*2^(M-L),做了一个换算下标都是N (0,0,0,0)(0,2,0,2)(0,1,2,3) */

/* 计算和转换因子乘积 */
TempReal2 = FFT_Input[ButterflyIndex2*2] * costab[ P ] +  FFT_Input[ButterflyIndex2*2+1] * sintab[ P ];
TempImag2 = FFT_Input[ButterflyIndex2*2+1] * costab[ P ] -  FFT_Input[ButterflyIndex2*2] * sintab[ P ] ;
TempReal1 = FFT_Input[ButterflyIndex1*2];
TempImag1 = FFT_Input[ButterflyIndex1*2+1];

/* 蝶形运算 */
FFT_Input[ButterflyIndex1*2] = TempReal1 + TempReal2;
FFT_Input[ButterflyIndex1*2+1] = TempImag1 + TempImag2;
FFT_Input[ButterflyIndex2*2] = TempReal1 - TempReal2;
FFT_Input[ButterflyIndex2*2+1] = TempImag1 - TempImag2;
}
}

Butterfly_NoPerGroup<<=1; /* 一组中蝶形的个数(1,2,4) */
}
}

/*
*********************************************************************************************************
* 函数名: Infinite_FFT
* 功能说明:
*    此函数用于对输入的复数组进行快速傅里叶变换（FFT），并完成后续的相关操作，可直接调用该函数。
*    不过需要注意的是，该函数处理后基波衰减相较于DSP库处理后的要大。
*    在使用此函数前，相关变量需在主函数中定义，示例如下：
*    Unlimited_FFT infinite_fft_s[x]
*    Unlimited_FFT infinite_fft_windowed_s[x]
*
* 形参:
*    - ADC_DMA_floatData: 指向需要进行FFT的数据的指针，数据类型为float32_t。
* - FFT_Input 指向用于存储FFT输入复数数组的指针，实部和虚部交替存储
* - FFT_Output 指向用于存储FFT输出幅频数组的指针
*    - Infinite_FFT_Length: 进行FFT运算的点数，数据类型为uint32_t。
*
* 使用举例:
* // FFT测试
* Infinite_FFT(ADC_Channels[0].DMA_floatData, infinite_fft_s[0].fft_Input, infinite_fft_s[0].fft_Output, MAX_FFT_N);
* Infinite_FFT(ADC_Channels[1].DMA_floatData, infinite_fft_s[1].fft_Input, infinite_fft_s[1].fft_Output, MAX_FFT_N);
*
* // 串口打印
* for(uint32_t i=0;i<ADC_DMA_SAMPLE_LENGTH;i++){
* printf("%f,%f,%f,%f\r\n", ADC_Channels[0].DMA_floatData*ZOOM, ADC_Channels[1].DMA_floatData*ZOOM, infinite_fft_s[0].fft_Output*4, infinite_fft_s[1].fft_Output*4);
* }
* // LCD
* lcd_clear(BLACK);
* lcd_draw_line(288, 0, 288, 800, RED);
* lcd_draw_line(96 , 0, 96 , 800, RED);
* for(uint32_t i=0;i<800;i++){
* lcd_draw_line(ADC_Channels[0].DMA_floatData*96.0/65535.0+384-48, i, ADC_Channels[0].DMA_floatData[i+1]*96.0/65535.0+384-48, i+1, YELLOW);
* lcd_draw_line(ADC_Channels[1].DMA_floatData*96.0/65535.0+192-48, i, ADC_Channels[1].DMA_floatData[i+1]*96.0/65535.0+192-48, i+1, LIGHTBLUE);
*
* lcd_draw_line(infinite_fft_s[0].fft_Output*50.0+288, i+10, infinite_fft_s[0].fft_Output[i+1]*50.0+288, i+10+1, RED);
* lcd_draw_line(infinite_fft_s[1].fft_Output*50.0+96, i+10, infinite_fft_s[1].fft_Output[i+1]*50.0+96, i+10+1, RED);
* }
*
* 返回值:
*    Infinite_FFT_End_Flag: FFT运算结束标志位，数据类型为uint8_t。当FFT运算完成时，该标志位被置为1并返回。
*
* 详细步骤:
*    1. 若MAX_FFT_N不为8192且不为16384，则调用InitTableFFT函数计算一批sin和cos系数。
*    2. 将输入的实数数据乘以缩放因子ZOOM后赋值给复数结构体的实部，虚部置为0。
*    3. 调用cfft函数进行指定点数的快速FFT运算。
*    4. 计算每个复数的幅频，即实部和虚部的平方和的平方根。
*    5. 对幅度谱进行归一化处理，直流分量除以FFT点数，其他频率分量除以FFT点数的一半。
**********************************************************************************************************
*/
uint8_t Infinite_FFT(float32_t *ADC_DMA_floatData, float32_t *FFT_Input, float32_t *FFT_Output, uint32_t Infinite_FFT_Length)
{
uint8_t Infinite_FFT_End_Flag = 0;
Infinite_FFT_End_Flag = 1;
/* 计算一批sin，cos系数 */
#if (MAX_FFT_N != 8192) && (MAX_FFT_N != 16384) && (MAX_FFT_N != 32768)
InitTableFFT(Infinite_FFT_Length);
#endif

for(uint32_t i=0; i < Infinite_FFT_Length; i++)
{
FFT_Input[i*2] = ADC_DMA_floatData * ZOOM;
FFT_Input[i*2+1] = 0;
}

/* MAX_FFT_N点快速FFT */
cfft(FFT_Input, Infinite_FFT_Length);

/* 计算幅频 */
// for(uint32_t i=0; i<Infinite_FFT_Length; i++)
// {
// FFT_Output = arm_sqrt_f32((float32_t)(FFT_Input[i*2] *FFT_Input[i*2]+ FFT_Input[i*2+1]*FFT_Input[i*2+1]), &FFT_Input[i*2]);
// }
arm_cmplx_mag_f32(FFT_Input, FFT_Output, Infinite_FFT_Length);

/* 对幅度谱进行归一化处理 */
/* 直流分量进行归一化 */
FFT_Output[0] /= Infinite_FFT_Length;
/* 其他频率分量进行归一化 */
for (uint32_t i = 1; i < Infinite_FFT_Length; i++) {
FFT_Output /= (Infinite_FFT_Length / 2.0);
}

return Infinite_FFT_End_Flag;
}

/**
* @brief 在FFT输出结果中查找基波幅值、索引、频率
*
* 该函数在FFT输出缓冲区中查找除直流分量外的最大幅值（峰值），记录其幅值和索引，
* 并根据采样率和FFT点数计算出对应的频率。
*
* @param FFT_OutputBuf 指向FFT输出结果的浮点数组指针，索引0为直流分量，函数从索引1开始查找
* @param fft_PeakAmp 指向用于存储FFT输出结果中峰值幅值的浮点变量指针
* @param fft_PeakAmpIndex 指向用于存储FFT输出结果中峰值幅值所在索引的无符号32位整数变量指针
* @param SP_rate 采样率，用于计算峰值对应的实际频率
* @param fft_n FFT点数，用于计算峰值对应的实际频率
*
* @return 峰值对应的实际频率
*/
float32_t FindfftPeak_f32(float32_t *FFT_OutputBuf, float32_t *fft_PeakAmp, uint32_t *fft_PeakAmpIndex, float32_t SP_rate, uint32_t fft_n)
{
float32_t f = 0;
arm_max_f32(&FFT_OutputBuf[1], fft_n/2, fft_PeakAmp, fft_PeakAmpIndex);
*fft_PeakAmpIndex += 1;
f = *fft_PeakAmpIndex * (SP_rate / fft_n);
return f;
}

/*
* 函数名: Find_PhaseAngle
* 功能说明: 求出该数组中幅值最大的位置对应的相位角
* 形参: ARR  ADC采样数组
*       N 遍历的总数量的一半
* 返回值: 相位角，单位：度
*/
float32_t Find_PhaseAngle(float32_t *FFT_Input, float32_t *FFT_OutputBuf, uint32_t fft_n)
{
float32_t pResult;
uint32_t pIndex;
float32_t phase_TEMP;
arm_max_f32(&FFT_OutputBuf[1], fft_n/2, &pResult, &pIndex);;//找出幅值的最大值下标
phase_TEMP = atan2f(FFT_Input[2 * pIndex+1], FFT_Input[2 * pIndex]) * 180.0f / 3.1415926f;//计算相位角
return phase_TEMP;
}

/**
* @brief 计算并去除输入信号直流分量
*
* 此函数计算 data 数组中信号的直流分量（平均值），并将其从信号中移除，
* 计算过程和去除操作中会使用 ZOOM 和 IZOOM 进行缩放。
*
* @param inputdata 指向存储信号数据的数组指针，函数会修改该数组内容
* @param DCpart 指向用于存储计算出的直流分量的变量指针
* @param length 数组 inputdata 的长度，即信号采样点数
*
* @return 无
*/
void RemoveDCOffset(float32_t* inputdata, float32_t* DCpart, uint32_t length)
{
float32_t sum = 0;
/* 计算信号的直流分量（平均值）*/
for (uint32_t i = 0; i < length; i++) sum += inputdata * ZOOM;
/* 计算平均值，即直流分量 */
*DCpart = sum / length;
/* 去除直流分量并进行缩放 */
for (uint32_t i = 0; i < length; i++) inputdata = (inputdata * ZOOM - *DCpart) * IZOOM;
}

sssuuu · 发表于 2025-3-19 18:06:34

最新测试：
从625Hz - 7500Hz，步进为625Hz来测试，最终计算得都是正确的频率
直到测8125Hz 算得 7500Hz
测8750Hz 算得8125Hz
测9375Hz 算得8750Hz
测10kHz 算得9375Hz
测10625Hz 算得10Hz
为什么会这样！
大佬们啊，救救我吧！

ME_Engineer · 发表于 2025-3-19 18:48:59

你分辨率搞错了吧？采样率640ksps=640000Hz，采样1024个点，频谱分辨率是640000/1024=625Hz。

sssuuu · 发表于 2025-3-20 11:16:03

ME_Engineer 发表于 2025-3-19 18:48
你分辨率搞错了吧？采样率640ksps=640000Hz，采样1024个点，频谱分辨率是640000/1024=625Hz。

谢谢大佬提醒，确实是靠定时器触发的ADC采样率是设置不准的，设置的640ksps，实际上是641ksps，然后就有很严重的频谱泄露，现在改用EXTI外部触发了，谢谢哥

风过不留痕 · 发表于 2025-3-20 13:31:40

sssuuu 发表于 2025-3-20 11:16
谢谢大佬提醒，确实是靠定时器触发的ADC采样率是设置不准的，设置的640ksps，实际上是641ksps，然后就有 ...

是改成外部给pwm信号了？信号发生器吗

sssuuu · 发表于 2025-3-20 14:23:50

风过不留痕发表于 2025-3-20 13:31
是改成外部给pwm信号了？信号发生器吗

对，用信号发生器给的方波信号，设置的采样率就很准

		自动登录	找回密码
密码			立即注册

[DSP] fft测量频率不准问题

浏览过的版块