搜索
bottom↓
回复: 187

使用STM32 的DSP库进行FFT变换说明及例程

  [复制链接]

出0入0汤圆

发表于 2010-8-11 13:40:09 | 显示全部楼层 |阅读模式
/*
*********************************************************************************************************
FileName:dsp_asm.h
*********************************************************************************************************
*/

#ifndef  __DSP_ASM_H__
#define  __DSP_ASM_H__
*********************************************************************************************************
*                                            FUNCTION PROTOTYPES
*********************************************************************************************************
*/

void    dsp_asm_test(void);
void  dsp_asm_init(void);

#endif                                                          /* End of module include.                               */
/*8888888888888888888888888888888888888888888888888888888888888888*/
/*8888888888888888888888888888888888888888888888888888888888888888*/
/*
* FileName:dsp_asm.c
* Author:Bobby.Chen
* Email:heroxx@163.com
* Date:2010-08-11
* Description:This file showes how to use the dsp library in mdk project.
*                  使用三角函数生成采样点,供FFT计算
*                  进行FFT测试时,按下面顺序调用函数即可:
*                   dsp_asm_init();
*                   dsp_asm_test();
*/
#include "stm32f10x.h"
#include "dsp_asm.h"
#include "stm32_dsp.h"
#include "table_fft.h"
#include <stdio.h>
#include <math.h>


/*
*********************************************************************************************************
*                                           LOCAL CONSTANTS
*********************************************************************************************************
*/
#define PI2 6.28318530717959
// Comment the lines that you don't want to use.
// 要模拟FFT,请注释掉其他的预定义
// 此处也可以全部注释掉,在MDK的工程属性->"C/C++"->"Preprocessor Symbols"-"Define:"中添加NPT_XXX项目
// 但是这样做法的缺点是每次修改XXX数据,都会导致MDK下次编译时会编译全部文件,速度太慢。
//#define NPT_64 64
#define NPT_256 256
//#define NPT_1024 1024

// N=64,Fs/N=50Hz,Max(Valid)=1600Hz
// 64点FFt,采样率3200Hz,频率分辨率50Hz,测量最大有效频率1600Hz
#ifdef NPT_64
#define NPT 64
#define Fs  3200
#endif

// N=256,Fs/N=25Hz,Max(Valid)=3200Hz
// 256点FFt,采样率6400Hz,频率分辨率25Hz,测量最大有效频率3200Hz
#ifdef NPT_256
#define NPT 256
#define Fs  6400
#endif

// N=1024,Fs/N=5Hz,Max(Valid)=2560Hz
// 1024点FFt,采样率5120Hz,频率分辨率5Hz,测量最大有效频率2560Hz
#ifdef NPT_1024
#define NPT 1024
#define Fs  5120
#endif

/*
*********************************************************************************************************
*                                       LOCAL GLOBAL VARIABLES
*********************************************************************************************************
*/
extern  uint16_t  TableFFT[];
long lBUFIN[NPT];         /* Complex input vector */
long lBUFOUT[NPT];        /* Complex output vector */
long lBUFMAG[NPT];/* Magnitude vector */
/*
*********************************************************************************************************
*                                      LOCAL FUNCTION PROTOTYPES
*********************************************************************************************************
*/
void dsp_asm_powerMag(void);

/*
*********************************************************************************************************
*   Initialize data tables for lBUFIN
* 模拟采样数据,采样数据中包含3种频率正弦波:50Hz,2500Hz,2550Hz
* lBUFIN数组中,每个单元数据高字(高16位)中存储采样数据的实部,低字(低16位)存储采样数据的虚部(总是为0)
*********************************************************************************************************
*/
void  dsp_asm_init()
{
  u16 i=0;
  float fx;
  for(i=0;i<NPT;i++)
  {
    fx  = 4000 * sin(PI2*i*50.0/Fs) + 4000 * sin(PI2*i*2500.0/Fs) + 4000*sin(PI2*i*2550.0/Fs);
    lBUFIN = ((s16)fx)<<16;
  }
}

/*
*********************************************************************************************************
*   Test FFT,calculate powermag
*   进行FFT变换,并计算各次谐波幅值
*********************************************************************************************************
*/
void  dsp_asm_test()
{
// 根据预定义选择合适的FFT函数
#ifdef NPT_64
  cr4_fft_64_stm32(lBUFOUT, lBUFIN, NPT);
#endif

#ifdef NPT_256
  cr4_fft_256_stm32(lBUFOUT, lBUFIN, NPT);
#endif

#ifdef NPT_1024
  cr4_fft_1024_stm32(lBUFOUT, lBUFIN, NPT);
#endif

  // 计算幅值
  dsp_asm_powerMag();

//  printf("No. Freq      Power\n");
//  for(i=0;i<NPT/2;i++)
//  {
//    printf("%4d,%4d,%10d,%10d,%10d\n",i,(u16)((float)i*Fs/NPT),lBUFMAG,(lBUFOUT>>16),(lBUFOUT&0xffff));
//  }
//  printf("*********END**********\r\n");
}
/*
*********************************************************************************************************
*   Calculate powermag
*   计算各次谐波幅值
*   先将lBUFOUT分解成实部(X)和虚部(Y),然后计算赋值(sqrt(X*X+Y*Y)
*********************************************************************************************************
*/
void dsp_asm_powerMag(void)
{
  s16 lX,lY;
  u32 i;
  for(i=0;i<NPT/2;i++)
  {
    lX  = (lBUFOUT << 16) >> 16;
    lY  = (lBUFOUT >> 16);
    {
    float X    = NPT * ((float)lX) /32768;
    float Y    = NPT * ((float)lY) /32768;
    float Mag = sqrt(X*X + Y*Y)/NPT;
    lBUFMAG    = (u32)(Mag * 65536);
    }
  }
}


// 笔者使用的是金牛开发板,CPU为STM32F107VC;JLink V8,MDK-ARM 4.10

// 注意FFT运算结果的对称性,也即256点的运算结果,只有前面128点的数据是有效可用的。
// 64点FFT运算结果图(局部):

FFT_64 (原文件名:1.JPG)

上图中,数组下标X对应的谐波频率为:N×Fs/64=N×3200/64=N*50Hz.

lBUFMAG[1] 对应 50Hz谐波幅值

上图中由于FFT分辨率50HZ,最大只能识别1600Hz谐波,导致结果中出现错误的数据。
// 256点FFT运算结果图(局部):


FFT_256 (原文件名:2.JPG)


上图中,数组下标X对应的谐波频率为:N×Fs/256=N×6400/256=N*25Hz.

lBUFMAG[2] 对应 2×25 =50Hz谐波幅值

lBUFMAG[100] 对应 100×25=2500Hz谐波幅值

lBUFMAG[102] 对应 102×25=2550Hz谐波幅值


// 1024点FFT运算结果图(局部):


FFT_1024 (原文件名:3.JPG)

上图中,数组下标X对应的谐波频率为:N×Fs/1024=N×5120/1024=N*5Hz.

lBUFMAG[10] 对应 10×5 =50Hz谐波幅值

lBUFMAG[500] 对应 500×5=2500Hz谐波幅值

lBUFMAG[510] 对应 510×5=2550Hz谐波幅值



FFT_Code示例代码FFT2.zip(文件大小:490K) (原文件名:FFT2.zip)

更正:
文件dsp_g2.c中dsp_asm_powerMag函数应按如下方式更改:
void dsp_asm_powerMag(void)
{
  s16 lX,lY;
  u32 i;
  for(i=0;i<NPT/2;i++)
  {
    lX  = (lBUFOUT << 16) >> 16;
    lY  = (lBUFOUT >> 16);
        {
        float X        = NPT * ((float)lX) /32768;
        float Y        = NPT * ((float)lY) /32768;
  float Mag = sqrt(X*X + Y*Y)/NPT;
  if(i==0)
    lBUFMAG        = (u32)(Mag * 32768);
  else
          lBUFMAG        = (u32)(Mag * 65536);
    }
  }
}
原因在于计算幅值时直流分量的幅值为sqrt(X*X+Y*Y)/NPT,不用除2

阿莫论坛20周年了!感谢大家的支持与爱护!!

如果想吃一顿饺子,就得从冰箱里取出肉,剁馅儿,倒面粉、揉面、醒面,擀成皮儿,下锅……
一整个繁琐流程,就是为了出锅时那一嘴滚烫流油的热饺子。

如果这个过程,禁不住饿,零食下肚了,饺子出锅时也就不香了……《非诚勿扰3》

出0入0汤圆

发表于 2010-8-11 13:45:35 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-8-11 14:23:02 | 显示全部楼层
有用

出0入42汤圆

发表于 2010-8-11 14:48:07 | 显示全部楼层
不错,mark

出0入0汤圆

发表于 2010-8-12 11:07:33 | 显示全部楼层
很好,适合初学者。

出0入0汤圆

发表于 2010-8-12 12:14:37 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-8-28 10:40:07 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-8-28 11:15:44 | 显示全部楼层
dsp_asm_powerMag()
这个函数实在废时间,比FFT运算时间多多多多

出0入0汤圆

 楼主| 发表于 2010-8-28 14:52:49 | 显示全部楼层
回复【7楼】wangguanfu  
dsp_asm_powermag()
这个函数实在废时间,比fft运算时间多多多多
-----------------------------------------------------------------------

确实如此,函数中应该采用定点数来计算,另外可以考虑用更优化的开放算法。

出0入0汤圆

发表于 2010-9-16 16:59:43 | 显示全部楼层
MARK

出0入0汤圆

发表于 2010-9-16 20:18:14 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-9-16 21:55:08 | 显示全部楼层
mark!~

出0入0汤圆

发表于 2010-9-16 22:31:17 | 显示全部楼层
mark???

出0入0汤圆

发表于 2010-9-16 22:56:53 | 显示全部楼层
很好谢谢

出0入0汤圆

发表于 2010-9-16 23:00:34 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-9-17 09:18:11 | 显示全部楼层
mark

出0入14汤圆

发表于 2010-9-17 09:35:17 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-9-17 09:41:28 | 显示全部楼层
jkkkkkkk

出0入0汤圆

发表于 2010-9-17 16:48:17 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-9-17 19:09:45 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-9-17 23:19:48 | 显示全部楼层
记号!

出0入0汤圆

发表于 2010-11-1 17:32:48 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-11-1 18:39:28 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-11-1 19:29:07 | 显示全部楼层
mark!~

出0入0汤圆

发表于 2010-11-1 20:33:00 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-11-23 15:21:58 | 显示全部楼层
不错,mark

出0入24汤圆

发表于 2010-11-24 15:29:02 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-11-24 21:11:06 | 显示全部楼层
学习

出0入0汤圆

发表于 2010-11-29 16:55:18 | 显示全部楼层
请教LZ 频率为0模值/NPT 是不是等于直流分量?

出0入0汤圆

发表于 2010-11-29 18:06:12 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-11-29 22:51:24 | 显示全部楼层
留下记号,这个是好东西呢。

出0入0汤圆

发表于 2010-11-30 00:07:20 | 显示全部楼层
正好让我赶上了,谢谢楼主!

出0入0汤圆

发表于 2010-12-2 15:51:44 | 显示全部楼层
记号!~这个要用到的....

出0入0汤圆

发表于 2010-12-3 10:39:16 | 显示全部楼层
mark!

出0入0汤圆

发表于 2010-12-3 15:28:28 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-12-3 16:59:38 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-12-11 21:57:17 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-12-11 22:12:04 | 显示全部楼层
mark

出0入0汤圆

发表于 2010-12-11 22:28:41 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-2-10 00:01:50 | 显示全部楼层
很好很好,学习一下,正用得上。

出0入0汤圆

发表于 2011-2-10 08:27:41 | 显示全部楼层
要是暑假看到就好了,fft想了就痛!

出0入0汤圆

发表于 2011-2-10 08:44:34 | 显示全部楼层
回复【40楼】lyzhangxiang
要是暑假看到就好了,fft想了就痛!
-----------------------------------------------------------------------

是的

出0入0汤圆

发表于 2011-3-12 21:37:29 | 显示全部楼层
请教楼主,我试过64点的可以,但是256的就不可以了,这是为什么呢

出0入0汤圆

发表于 2011-3-12 21:44:39 | 显示全部楼层
标记

出0入0汤圆

发表于 2011-3-13 18:21:34 | 显示全部楼层
学习了,好东西啊!

出0入0汤圆

发表于 2011-3-13 22:02:34 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-5-13 16:20:18 | 显示全部楼层
做1024点的,一直有问题。
还为找到原因?
FFT前的序列输入的幅值有要求吗?比如1000~2000?

出0入0汤圆

发表于 2011-6-14 16:37:26 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-6-22 20:39:59 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-7-16 20:13:07 | 显示全部楼层
回复【楼主位】moyuker  
-----------------------------------------------------------------------
MARK

出0入0汤圆

 楼主| 发表于 2011-7-21 11:52:39 | 显示全部楼层
程序更新了ADC采样部分。

程序源文件,请下载后将文件后缀改为7z。ourdev_659734CIMQQR.zip(文件大小:218K) (原文件名:FFT.zip)

主要代码如下:
/*        根据采样频率设定定时器        */
void    adc_TIM_Init(void)
{
        TIM_TimeBaseInitTypeDef   TIM_TimeBaseStructure;
//        TIM_OCInitTypeDef         TIM_OCInitStructure;

        RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM3,ENABLE);

        /* Time Base configuration */
        TIM_TimeBaseStructInit(&TIM_TimeBaseStructure);
#ifdef        NPT_64
        TIM_TimeBaseStructure.TIM_Prescaler = 0;      
        TIM_TimeBaseStructure.TIM_Period = 22499;         
#endif
#ifdef        NPT_256
        TIM_TimeBaseStructure.TIM_Prescaler = 0;      
        TIM_TimeBaseStructure.TIM_Period = 11249;         
#endif
#ifdef        NPT_1024
        TIM_TimeBaseStructure.TIM_Prescaler = 0;      
        TIM_TimeBaseStructure.TIM_Period = 14062;         
#endif
        TIM_TimeBaseStructure.TIM_ClockDivision = 0x0;   
        TIM_TimeBaseStructure.TIM_CounterMode = TIM_CounterMode_Up;  
        TIM_TimeBaseInit(TIM3, &TIM_TimeBaseStructure);

        TIM_SelectOutputTrigger(TIM3,TIM_TRGOSource_Update);
        TIM_Cmd(TIM3,ENABLE);
        TIM_ITConfig(TIM3,TIM_IT_Update,ENABLE);
}

void    ac_adc_Init(void)
{
        ADC_InitTypeDef ADC_InitStructure;
        DMA_InitTypeDef DMA_InitStructure;


        ADC_DeInit(ADC1);
        ADC_DeInit(ADC2);
        RCC_APB2PeriphClockCmd(RCC_APB2Periph_ADC1 | RCC_APB2Periph_ADC2,ENABLE);
        RCC_AHBPeriphClockCmd(RCC_AHBPeriph_DMA1,ENABLE);

        /* DMA1 channel1 configuration ----------------------------------------------*/
        DMA_DeInit(DMA1_Channel1);
        DMA_InitStructure.DMA_PeripheralBaseAddr = (uint32_t)ADC1_DR_Address;
        DMA_InitStructure.DMA_MemoryBaseAddr = (uint32_t)ADC_DualConvertedValueTab;
        DMA_InitStructure.DMA_DIR = DMA_DIR_PeripheralSRC;
        DMA_InitStructure.DMA_BufferSize = NPT;
        DMA_InitStructure.DMA_PeripheralInc = DMA_PeripheralInc_Disable;
        DMA_InitStructure.DMA_MemoryInc = DMA_MemoryInc_Enable;
        DMA_InitStructure.DMA_PeripheralDataSize = DMA_PeripheralDataSize_Word;
        DMA_InitStructure.DMA_MemoryDataSize = DMA_MemoryDataSize_Word;
        DMA_InitStructure.DMA_Mode = DMA_Mode_Circular;//DMA_Mode_Normal;
        DMA_InitStructure.DMA_Priority = DMA_Priority_High;
        DMA_InitStructure.DMA_M2M = DMA_M2M_Disable;
        DMA_Init(DMA1_Channel1, &DMA_InitStructure);

        DMA_ITConfig(DMA1_Channel1,DMA1_IT_TC1,ENABLE);
        DMA_Cmd(DMA1_Channel1, ENABLE);

        /* ADC1 configuration ------------------------------------------------------*/
        ADC_InitStructure.ADC_Mode  = ADC_Mode_RegSimult;
        ADC_InitStructure.ADC_ScanConvMode = DISABLE;
        ADC_InitStructure.ADC_ContinuousConvMode = DISABLE;
        ADC_InitStructure.ADC_ExternalTrigConv = ADC_ExternalTrigConv_T3_TRGO;//ADC_ExternalTrigConv_None;
        ADC_InitStructure.ADC_DataAlign = ADC_DataAlign_Right;
        ADC_InitStructure.ADC_NbrOfChannel = 1;
        ADC_Init(ADC1, &ADC_InitStructure);
        ADC_RegularChannelConfig(ADC1, ADC_Channel_0, 1, ADC_SampleTime_239Cycles5);   
        ADC_DMACmd(ADC1, ENABLE);

        /* ADC2 configuration ------------------------------------------------------*/
        ADC_InitStructure.ADC_Mode = ADC_Mode_RegSimult;
        ADC_InitStructure.ADC_ScanConvMode = DISABLE;
        ADC_InitStructure.ADC_ContinuousConvMode = DISABLE;
        ADC_InitStructure.ADC_ExternalTrigConv = ADC_ExternalTrigConv_T3_TRGO;//ADC_ExternalTrigConv_None;
        ADC_InitStructure.ADC_DataAlign = ADC_DataAlign_Right;
        ADC_InitStructure.ADC_NbrOfChannel = 1;
        ADC_Init(ADC2, &ADC_InitStructure);
        ADC_RegularChannelConfig(ADC2, ADC_Channel_1, 1, ADC_SampleTime_239Cycles5);
        ADC_ExternalTrigConvCmd(ADC2, ENABLE);


        ADC_Cmd(ADC1, ENABLE);
        ADC_ResetCalibration(ADC1);
        while (ADC_GetResetCalibrationStatus(ADC1));
        ADC_StartCalibration(ADC1);
        while (ADC_GetCalibrationStatus(ADC1));

        ADC_Cmd(ADC2, ENABLE);
        ADC_ResetCalibration(ADC2);
        while (ADC_GetResetCalibrationStatus(ADC2));
        ADC_StartCalibration(ADC2);
        while (ADC_GetCalibrationStatus(ADC2));

        ADC_ITConfig(ADC1,ADC_IT_EOC,ENABLE);

}


void    DMA1_Channel1_IRQHandler(void)
{
        if (DMA_GetITStatus(DMA1_IT_TC1) != RESET) {
                LED5_Switch();
                ui_data_to_dsp(ADC_DualConvertedValueTab);
                DMA_ClearITPendingBit(DMA1_IT_TC1);
        }

        if (DMA_GetITStatus(DMA1_IT_TE1) != RESET) {
                DMA_ClearITPendingBit(DMA1_IT_TE1);
        }
}

void        ui_data_to_dsp(uint32_t* ADC_VTab)
{
        if(_bFlagNeedDSP != 0) {
                return;
        }
#ifndef DSP_TEST
        memcpy(ADC_ValueTable,ADC_VTab,sizeof(uint32_t) * NPT);
#endif

        _bFlagNeedDSP        = 1;
}

出0入0汤圆

发表于 2011-7-21 12:50:56 | 显示全部楼层
mark 参考一下

出0入0汤圆

发表于 2011-7-21 13:12:29 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-7-21 13:46:56 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-7-21 21:11:24 | 显示全部楼层
以后应该会用到……

出0入8汤圆

发表于 2011-7-21 21:31:42 | 显示全部楼层
good

出0入0汤圆

发表于 2011-7-21 22:11:26 | 显示全部楼层
标记一下

出0入0汤圆

发表于 2011-7-22 07:54:29 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-7-22 11:25:29 | 显示全部楼层
mark!!!!

出0入0汤圆

发表于 2011-7-31 20:04:14 | 显示全部楼层
回复【50楼】moyuker
-----------------------------------------------------------------------

你好,最近小弟看了下你的程序,有个问题,是不是stm32库做fft处理受频率限制,最高频率也是用256点时的3200hz啊??

出0入0汤圆

发表于 2011-7-31 20:40:33 | 显示全部楼层
没用过stm的库,请教一下:输入数据加窗吗?汉明窗还是其他什么?

出0入0汤圆

 楼主| 发表于 2011-8-1 23:15:14 | 显示全部楼层
回复【60楼】bingguji  
回复【50楼】moyuker
-----------------------------------------------------------------------
你好,最近小弟看了下你的程序,有个问题,是不是stm32库做fft处理受频率限制,最高频率也是用256点时的3200hz啊??
-----------------------------------------------------------------------

频率分辨率=采样频率/FFT点数
256点FFt,采样率6400Hz
频率分辨率6400/256=25Hz,测量最大有效频率6400/2=3200Hz

出0入0汤圆

 楼主| 发表于 2011-8-1 23:16:03 | 显示全部楼层
回复【61楼】leoyang  
没用过stm的库,请教一下:输入数据加窗吗?汉明窗还是其他什么?
-----------------------------------------------------------------------

此程序中没有进行加窗处理。

出0入0汤圆

发表于 2011-8-2 20:51:09 | 显示全部楼层
下载了LZ的程序后,发现里面还有个LZ自己也的FFT函数啊,转换的精度怎么样呢?
在测试程序里面
用DSP库的那个程序
void  dsp_asm_init()
{
  u16 i=0;
  float fx;
  for(i=0;i<NPT;i++)
  {
    fx  = 4000 * sin(PI2*i*50.0/Fs) + 4000 * sin(PI2*i*2500.0/Fs) + 4000*sin(PI2*i*2550.0/Fs);
    lBUFIN = ((s16)fx)<<16;
  }
}

你把这个改成下面这这个函数
void  dsp_asm_init()
{
  u16 i=0;
  float fx;
  u16 NS;
  NS=0;
  for(i=0;i<NPT;i++)
  {
    fx  = 4000 * sin(PI2*i*50.0/Fs+NS*PI2/12);
    lBUFIN = ((s16)fx)<<16;
  }
}
就是只给采样函数加了初相,在NS取值不同的时候(我测试的是NS=0~12,刚好360度)是后,进行FFT的输出的值是不同的

出0入0汤圆

 楼主| 发表于 2011-8-2 22:18:10 | 显示全部楼层
回复【64楼】dwlovework  
-----------------------------------------------------------------------
C语言的基2 FFT函数是从网上Copy下来的,准确率很高,但是计算速度慢。
我测试的结果
        
点击图片看大图
(原文件名:未标题-1.jpg)


这个现在改成STM32F103了。ourdev_663680IUER4E.zip(文件大小:564K) (原文件名:FFT2.zip)

出0入0汤圆

发表于 2011-8-2 23:56:23 | 显示全部楼层
FFT 应该是DSP的专长,用MCU去实现,运算速度不慢才怪。

出15入9汤圆

发表于 2011-8-3 00:51:37 | 显示全部楼层
这个要关注,呵呵。

出0入0汤圆

发表于 2011-8-3 10:07:57 | 显示全部楼层
回复【65楼】moyuker
-----------------------------------------------------------------------
那LZ有没有测试如果给信号加上初相,每一次采样的初相都不相同的话,使用DSP库的64点fft计算结果会不一样把,其余不变的情况下,初相改变,求取的幅值最大差16。。。。这个问题怎么解决啊,加窗么?

出0入0汤圆

 楼主| 发表于 2011-8-3 10:13:20 | 显示全部楼层
回复【68楼】dwlovework  
-----------------------------------------------------------------------
main()
{
.....
  while(1)
  {

        dsp_asm_init();
          dsp_asm_test();
  }
}

u16 NS        = 0;
void  dsp_asm_init()
{
  u16 i=0;
  float fx;
  
  NS        = (NS >=12 ) ? 0 : NS;
  
  for(i=0;i<NPT;i++)
  {
    //fx  = 4000 * sin(PI2*i*50.0/Fs) + 4000 * sin(PI2*i*2500.0/Fs) + 4000*sin(PI2*i*2550.0/Fs);
    fx        = 4000 * sin(PI2*i*50.0/Fs+NS*PI2/12);
    lBUFIN = ((s16)fx)<<16;
  }

  NS++;
}

出0入0汤圆

发表于 2011-8-3 10:17:06 | 显示全部楼层
回复【65楼】moyuker
-----------------------------------------------------------------------

基2的那个FFT的算法也有这种问题么?初相改变,然后求取的幅值会波动。现在想给他做校正,也不知道该怎么着手才好,都搞了4.5天了,一点进展都没有。LZ有想法么?可不可以指点下

出0入0汤圆

 楼主| 发表于 2011-8-3 10:42:01 | 显示全部楼层
回复【70楼】dwlovework  
-----------------------------------------------------------------------

现在只采样了20ms一个周期的数据做64点FFT,效果会比较差。
你可以采样40ms做一次256点的FFT,然后将200ms内的连续5次FFT结果做平均。

至于加窗算法,之前只是实验了一下汉宁窗,但是没有进行实际的测量计算,不知道效果如何。

百度文库可以搜到很多这方面的论文。

http://wenku.baidu.com/view/5919d4fbaef8941ea76e0524.html

出0入0汤圆

发表于 2011-8-3 11:00:50 | 显示全部楼层
回复【65楼】moyuker
-----------------------------------------------------------------------

刚下了LZ的程序,用基2的FFT测试了下,结果是很准确,几乎不受NS取值的影响。计算时间是将近6ms,感觉还行

那么官方的DSP库中的FFT为什么会有这个问题的存在呢?不解。。。

出0入0汤圆

发表于 2011-8-5 15:02:00 | 显示全部楼层
细细看下

出0入0汤圆

发表于 2011-8-6 16:54:41 | 显示全部楼层
仔细看看

出0入0汤圆

发表于 2011-8-9 12:19:19 | 显示全部楼层
回复【楼主位】moyuker
-----------------------------------------------------------------------

楼主好,解说很详细,超赞。
最新的DSP库(UM0585 Rev2)出来了:
因为在旧版本的库上有个问题,不知楼主注意没有,就是例程跟说明文档对于实部,虚部位置的解说不一致。
新版本已经一致化了:
“All the signal samples must be 32-bit data containing the 16-bit real part followed by the
16-bit imaginary part (in the little Endian order: imaginary_real)”
也就是
低16bit 实部
高16bit 虚部

(这样一来才比较合理,X、cos……一般代表实部,Y、sin……一般代表虚部)

楼主加油出好帖哇哈,因为你的帖子帮到了一批人哇

出0入0汤圆

发表于 2011-8-10 11:45:59 | 显示全部楼层
真是好东西 方便开发 支持一个

出0入0汤圆

发表于 2011-8-11 02:39:25 | 显示全部楼层
回复【62楼】moyuker
-----------------------------------------------------------------------

那么那个采样频率是固定的么????

出0入0汤圆

发表于 2011-8-11 09:08:32 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-8-12 23:17:42 | 显示全部楼层
mark!!

出0入0汤圆

发表于 2011-8-14 16:01:31 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-8-22 10:33:27 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-8-22 22:24:23 | 显示全部楼层
我想问下楼主,得出了各次谐波幅值后,信号的幅值取基波的幅值还是各次谐波的平方和开根号??

出0入0汤圆

发表于 2011-8-23 13:33:45 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-8-23 17:40:39 | 显示全部楼层
回复【85楼】chwei
-----------------------------------------------------------------------

可想学,但是FFT原理不懂

出0入0汤圆

发表于 2011-8-24 22:06:32 | 显示全部楼层
怎么没人回答呢?

出0入0汤圆

发表于 2011-9-16 22:08:17 | 显示全部楼层
Mark。。。占座。。值得参考

出0入0汤圆

发表于 2011-10-4 14:22:52 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-12-11 21:30:55 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-12-11 21:39:38 | 显示全部楼层
有没人测试一下楼主写的这些,跟ARM出的CMSIS库里面的dsplib的代码效率pk结果?

出0入0汤圆

发表于 2011-12-18 11:49:30 | 显示全部楼层
还是不理解模值运算那部分,不就是一个复数求模吗,实部平方加虚部平方再开方不就完了,为什么搞得这么复杂,为什么直流分量还要除以2,还有就是输入为什么还是带虚部的,实际中输入的不都是实数吗,多此一举啦

出0入0汤圆

发表于 2011-12-18 12:35:49 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-12-18 23:55:05 | 显示全部楼层
很好!学习下。

出0入0汤圆

发表于 2011-12-19 20:37:09 | 显示全部楼层
mark

出0入0汤圆

发表于 2011-12-20 13:36:59 | 显示全部楼层
恩,原来是FFT后,输出的虚数和实数都是带符号位的,这样求模值就没有我想象的那么简单了

出0入0汤圆

发表于 2012-2-1 13:25:45 | 显示全部楼层
mark

出0入0汤圆

发表于 2012-2-28 21:19:20 | 显示全部楼层
楼主:
不管信号是从哪个相位开始采集的,只要符合64点 256点、、、、它的频谱还是保持一样。这不正是从频域角度去分析信号嘛。我还沉迷于怎么从0相位开始触发ADC的 呵呵谢谢楼主开题,学了不少东西

出0入0汤圆

发表于 2012-2-28 21:50:56 | 显示全部楼层
请教moyuker

你说的/*        获取交流信号中心线电压码值        */是什么意思?我理解这是在获取一个相对稳定的起点电压值,
∵实际上是在  /*等待ADC1采样稳定*/
        while(i-->0) {
                _base        = ADC_GetConversionValue(ADC1);
        }

期待楼主回复、、、

出0入0汤圆

发表于 2012-3-29 10:19:49 | 显示全部楼层
标记一下

出0入0汤圆

发表于 2012-5-8 11:46:58 | 显示全部楼层
顶,有个问题请教一下楼主:
for(i=0;i<NPT/2;i++)
  {
    lX  = (lBUFOUT << 16) >> 16;
    lY  = (lBUFOUT >> 16);
    {
    float X    = NPT * ((float)lX) /32768;
    float Y    = NPT * ((float)lY) /32768;
    float Mag = sqrt(X*X + Y*Y)/NPT;
    lBUFMAG    = (u32)(Mag * 65536);
中cr4_fft_64_stm32(DSP_OUT,DSP_IN,64);函数输出的DSP_OUT哪个是实部哪个是虚部??
官方文档中有
All the signal samples must be 32-bit data containg the 16-bit real part followed by the
16-bit imaginary part (in the little Endian order: imaginary_real).
结合文档中下面的例程,不理解上面这话的意思,
你怎么说* lBUFIN数组中,每个单元数据高字(高16位)中存储采样数据的实部,低字(低16位)存储采样数据的虚部(总是为0)呢 求解?
另外楼主有没有做过--->由各次谐波信息计算各次谐波对应的相位?

出0入0汤圆

发表于 2012-5-9 16:21:31 | 显示全部楼层
太耗时间了吧
回帖提示: 反政府言论将被立即封锁ID 在按“提交”前,请自问一下:我这样表达会给举报吗,会给自己惹麻烦吗? 另外:尽量不要使用Mark、顶等没有意义的回复。不得大量使用大字体和彩色字。【本论坛不允许直接上传手机拍摄图片,浪费大家下载带宽和论坛服务器空间,请压缩后(图片小于1兆)才上传。压缩方法可以在微信里面发给自己(不要勾选“原图),然后下载,就能得到压缩后的图片】。另外,手机版只能上传图片,要上传附件需要切换到电脑版(不需要使用电脑,手机上切换到电脑版就行,页面底部)。
您需要登录后才可以回帖 登录 | 注册

本版积分规则

手机版|Archiver|amobbs.com 阿莫电子技术论坛 ( 粤ICP备2022115958号, 版权所有:东莞阿莫电子贸易商行 创办于2004年 (公安交互式论坛备案:44190002001997 ) )

GMT+8, 2024-3-28 19:44

© Since 2004 www.amobbs.com, 原www.ourdev.cn, 原www.ouravr.com

快速回复 返回顶部 返回列表