dsPIC33AKシリーズ⑰ DSPの使い方2（FFT）

dsPIC33Aシリーズ
2024.08.26

概要

以前の記事にも有りますように、dsPICシリーズは「dsP」と頭文字がつくように、DSPが内蔵されています。今回は、dsPIC33AKシリーズのDSPを用いてFFT(*1)演算を実行し、従来との性能差を検証します。

FFTとは

FFTとは高速フーリエ変換(Fast Fourier Transform)の略。フーリエ変換とは時間領域の関数を周波数領域の関数に変換することを指し、それを高速に演算するアルゴリズムの事

関連記事

記事リンク
第1回dsPIC33AKシリーズについてdsPIC33Aシリーズに関して – ぴくおの電子工作的な何かWP (electricpico.com)
第2回開発ボード、コンフィグレーション設定、クロック設定についてdsPIC33AKシリーズ②開発ボード,コンフィグレーション設定,クロック設定について – ぴくおの電子工作的な何かWP (electricpico.com)
第3回CPU性能についてdsPIC33AKシリーズ③CPU性能について – ぴくおの電子工作的な何かWP (electricpico.com)
第4回FPU性能についてdsPIC33AKシリーズ④FPU性能について – ぴくおの電子工作的な何かWP (electricpico.com)
第5回DSPについてdsPIC33AKシリーズ⑤DSPについて – ぴくおの電子工作的な何かWP (electricpico.com)
第6回タイマー1割り込みの使い方dsPIC33AKシリーズ⑥Timer1割り込みの使い方 – ぴくおの電子工作的な何かWP (electricpico.com)
第7回PMUの使い方dsPIC33AKシリーズ⑦PMUの使い方 – ぴくおの電子工作的な何かWP (electricpico.com)
第8回高速ADCの使い方dsPIC33AKシリーズ⑧高速ADCの使い方 – ぴくおの電子工作的な何かWP (electricpico.com)
第9回高速OPAMPの使い方dsPIC33AKシリーズ⑨高速OPAMPの使い方 – ぴくおの電子工作的な何かWP (electricpico.com)
第10回DMAモジュールの使い方dsPIC33AKシリーズ⑩DMAモジュールの使い方 – ぴくおの電子工作的な何かWP (electricpico.com)
第11回WDTモジュールの使い方dsPIC33AKシリーズ⑪WDTモジュールの使い方 – ぴくおの電子工作的な何かWP (electricpico.com)
第12回DMTモジュールの使い方dsPIC33AKシリーズ⑫DMTモジュールの使い方 – ぴくおの電子工作的な何かWP (electricpico.com)
第13回I/Oインテグリティモニタモジュールの使い方dsPIC33AKシリーズ⑬I/Oインテグリティモニターモジュールの使い方 – ぴくおの電子工作的な何かWP (electricpico.com)
第14回QEIモジュールの使い方dsPIC33AKシリーズ⑭QEIモジュールの使い方 – ぴくおの電子工作的な何かWP (electricpico.com)
第15回UARTモジュールの使い方dsPIC33AKシリーズ⑮UARTモジュールの使い方 – ぴくおの電子工作的な何かWP (electricpico.com)
第16回DSPの使い方1 （デジタルフィルタ）dsPIC33AKシリーズ⑯ DSPの使い方1（デジタルフィルタ） – ぴくおの電子工作的な何かWP (electricpico.com)
第17回（本記事）DSPの使い方2 (FFT)dsPIC33AKシリーズ⑰ DSPの使い方2（FFT） – ぴくおの電子工作的な何かWP (electricpico.com)
(fig.1-1)各関連記事リンク
関連記事リンク
DSPライブラリの使い方1（デジタルフィルタ）DSPライブラリの使い方1（デジタルフィルタ） – ぴくおの電子工作的な何かWP (electricpico.com)
DSPライブラリの使い方2（FFT）DSPライブラリの使い方2（FFT） – ぴくおの電子工作的な何かWP (electricpico.com)
(fig.1-2)各関連記事リンク

開発環境

開発環境を以下に示します。
今回はCPUボード(EV02G02A)で動作確認を行います。

項目値リンク
ベースボードdsPIC33A CURIOSITY PLATFORM
DEVELOPMENT BOARD		dsPIC33A Curiosity Platform Development Board User’s Guide (microchip.com)
CPUボード(EV68M17A)EV68M17A – dsPIC33AK128MC106 Motor Control DIMdsPIC33AK128MC106 Motor Control Dual In-Line Module (DIM) Information Sheet (microchip.com)
CPUボード(EV02G02A)dsPIC33AK128MC106 General Purpose Dual In-Line Module (DIM)dsPIC33AK128MC106 General Purpose Dual In-Line Module (DIM) | Microchip Technology
統合開発環境MPLAB X IDE v6.20MPLAB® X IDE | Microchip Technology
コンパイラMPLAB XC DSC v3.10MPLAB® XC DSC Compiler | Microchip Technology
(fig.2-1)本記事の動作確認環境
(fig.2-2)動作環境

動作確認

ソースコード

インクルードファイル

コンフィグレーションファイル、クロック設定ファイルは以下のファイルをインクルードしてください。

■コンフィグレーションファイル　(config.h)

config.hダウンロード

■クロック設定ソースファイル(Clock_Driver.c)

Clock_Driver.cダウンロード

■クロックヘッダーファイル(clock_driver.c)

clock_driver.hダウンロード

ライブラリのインポート

DSPを使用する場合、dsp.hのインクルードとlibdsp-elf.aのインポートが必須です。

ソースコード全体

ダミーデータとして1kHzのサイン波を16kHzでサンプリングした（と仮定した）128ポイントのデータをfr_InputDataに格納します。

// 入力バッファ　128サンプル 1000Hz/16kHz fs
fractional	fr_InputData[ FFT_POINT ] = 
{
0	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
0	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207	,
-1	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
0	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207	,
-1	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
0	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207	,
-1	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
0	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207	,
-1	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
0	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207	,
-1	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
-1	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207	,
-1	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
-1	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207	,
-1	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
-1	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207

};

初期設定として以下のコードのように、回転因子の初期化と窓関数の初期化を行います。

/*------------------------------------------------------------------------*/
//回転因子の初期化
/*------------------------------------------------------------------------*/
　　tFactor = TwidFactorInit(
	FFT_POWEROF2,
	(fractcomplex *)twiddleFactors,
	0);
/*------------------------------------------------------------------------*/
//窓関数の初期化
/*------------------------------------------------------------------------*/
　　HanningInit(FFT_POINT,fr_Window);

初期化が終わった時点でFFTを実行します。

Calc_FFT();

■ソースコード全体

#include <stdio.h>
#include <stdlib.h>
#include <xc.h>
#include "config.h"
#include "clock_driver.h"
#include <dsp.h>
/*----------------------------------------------------------------------------*/
/* 定数定義*/
/*----------------------------------------------------------------------------*/
#define FFT_POWEROF2    7u
#define FFT_POINT       ( 1 << FFT_POWEROF2)
/*----------------------------------------------------------------------------*/
/*変数定義*/
/*----------------------------------------------------------------------------*/
//回転因子
fractcomplex twiddleFactors[FFT_POINT] __attribute__ ((section (".ybss, bss, ymemory"), aligned (FFT_POINT)));
// スペクトル(複素数) 
fractcomplex frc_FftData[ FFT_POINT ] __attribute__ ((section (".xbss, bss, xmemory"), aligned (FFT_POINT)));

//回転因子ポインタ
fractcomplex *tFactor = (fractcomplex * )&twiddleFactors;

// 入力バッファ　128サンプル 1000Hz/16kHz fs
fractional	fr_InputData[ FFT_POINT ] = 
{
0	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
0	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207	,
-1	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
0	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207	,
-1	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
0	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207	,
-1	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
0	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207	,
-1	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
0	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207	,
-1	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
-1	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207	,
-1	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
-1	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207	,
-1	,
410903206	,
759250124	,
992008093	,
1073741823	,
992008093	,
759250124	,
410903206	,
-1	,
-410903207	,
-759250125	,
-992008094	,
-1073741824	,
-992008094	,
-759250125	,
-410903207

};

// パワースペクトル（実数）
fractional	fr_PowerSpec[ FFT_POINT/2 ];	
fractional	fr_PowerSpec1[ FFT_POINT/2 ];
// 窓関数用
fractional fr_Window[ FFT_POINT ];
fractional fr_AfterWindow[ FFT_POINT ];
int32_t i4g_Time[16];


void Calc_FFT(void);
static void Uart_Send_Blocking_TX_Data(uint8_t *u1pa_Data);

/*----------------------------------------------------------------------------*/
/**
* @fn          main(int argc, char** argv) 
* @brief       PROJECT_4_1_1_DSP_FIR
* @param[in]   argc argument count
* @param[in]   argv argument vector
* @retval      EXIT_SUCCESS 成功
* @retval      EXIT_FAILURE 失敗
* @detail      DSPによるFFT
* @note        
 */
/*----------------------------------------------------------------------------*/
int main(int argc, char** argv) 
{
    uint8_t TxTemp[64];
    /*------------------------------------------------------------------------*/
    /*　クロック初期化*/
    /*------------------------------------------------------------------------*/
        vdg_Clock_Set_Register();    /* クロック初期化 */
  　/*-----------------------------------------------------------------------*/
    /* ピン設定*/
    /*-----------------------------------------------------------------------*/
	RPOR12bits.RP50R = 9u;		//UART1出力 RP50
	RPINR9bits.U1RXR = 52u;		//UART1入力 RP52
    /*-----------------------------------------------------------------------*/
    /* UART初期化*/
    /*-----------------------------------------------------------------------*/
	U1CON = 0x00000000u;
	U1STAT = 0x00000000u;
	U1BRG = (uint32_t) (100000000 / 115200);
	U1RXB = 0x00000000u;
	U1TXB = 0x00000000u;
	U1PA = 0x00000000u;
	U1PB = 0x00000000u;
	U1CHK = 0x00000000u;
	U1SCCON = 0x00000000u;
	U1UIR = 0x00000000u;

	U1CONbits.CLKMOD = 1u;
	U1CONbits.RXEN = 1u;	//受信許可
	U1CONbits.TXEN = 1u;	//送信許可
	U1CONbits.ON = 1u;
    /*-----------------------------------------------------------------------*/
    /* タイマー初期化*/
    /*-----------------------------------------------------------------------*/
	T1CONbits.ON = 1;	
	TMR1 = 0;
		
	i4g_Time[0] = TMR1;
  　/*------------------------------------------------------------------------*/
  　/* ①回転因子の初期化*/
  　/*------------------------------------------------------------------------*/
	tFactor = TwidFactorInit(
		FFT_POWEROF2,
		(fractcomplex *)twiddleFactors,
		0);
	i4g_Time[1] = TMR1;
  　/*------------------------------------------------------------------------*/
  　/* ②窓関数の初期化*/
  　/*------------------------------------------------------------------------*/
	HanningInit(FFT_POINT,fr_Window);
	i4g_Time[2] = TMR1;
  　/*------------------------------------------------------------------------*/
  　/* FFTの実行*/
  　/*------------------------------------------------------------------------*/
	Calc_FFT();
		
  　/*------------------------------------------------------------------------*/
  　/* 結果の出力*/
  　/*------------------------------------------------------------------------*/
	for (int i = 0;i< FFT_POINT;i ++)
	{
sprintf(TxTemp,"%d,%d,%d,%d\r\n",i,fr_InputData[i],fr_Window[i],fr_AfterWindow[i]);
		Uart_Send_Blocking_TX_Data(TxTemp);
	}

	for (int i = 0;i< FFT_POINT / 2 ;i ++)
	{
sprintf(TxTemp,"%d,%d,%d,%d,%d\r\n",i,frc_FftData[i].real,frc_FftData[i].imag,fr_PowerSpec[i],fr_PowerSpec1[i]);
		Uart_Send_Blocking_TX_Data(TxTemp);
	}

	while (1)
	{
         ;
	}
    return EXIT_SUCCESS;
}
/*----------------------------------------------------------------------------*/
/**
* @fn         Uart_Send_Blocking_TX_Data(uint8_t *u1pa_Data)
* @brief      UART送信の実行
* @param      *u1pa_Data 送信データポインタ
* @note        
*/
/*----------------------------------------------------------------------------*/
static void Uart_Send_Blocking_TX_Data(uint8_t *u1pa_Data)
{
　　while (*u1pa_Data != NULL)
　　{
	while(U1STATbits.TXBF == 1u){;}
	U1TXB = *u1pa_Data ++;
　　}	
}
/*----------------------------------------------------------------------------*/
/**
* @fn          Calc_FFT(void)
* @brief       FFTの実行
* @param　　　　なし
* @note        
*/
/*----------------------------------------------------------------------------*/
void Calc_FFT(void)
{
    int i;
    
  /*------------------------------------------------------------------------*/
  /* ③窓関数の適用（DSP関数）*/
  /*------------------------------------------------------------------------*/
	VectorWindow(FFT_POINT,
		(fractional *)fr_AfterWindow,
		(fractional *)fr_InputData,
		(fractional *)fr_Window);
	i4g_Time[3] = TMR1;
  /*------------------------------------------------------------------------*/
  /* ④入力データを複素数の実数に代入*/
  /*------------------------------------------------------------------------*/
	for (i = 0; i < FFT_POINT;i++)
	{
		frc_FftData[i].real = fr_AfterWindow[i];
		frc_FftData[i].imag = 0;
	}

	i4g_Time[4] = TMR1;
  /*------------------------------------------------------------------------*/
  /* ⑤FFTの実行（DSP関数）*/
  /*------------------------------------------------------------------------*/
	FFTComplexIP( 
		FFT_POWEROF2, 
		(fractcomplex *)frc_FftData, 
		(fractcomplex *)twiddleFactors,
		0);
	i4g_Time[5] = TMR1;
	
  /*------------------------------------------------------------------------*/
  /* ⑥ビット逆順ソート（DSP関数）*/
  /*------------------------------------------------------------------------*/
	BitReverseComplex( FFT_POWEROF2, frc_FftData ); 
	i4g_Time[6] = TMR1;
  /*------------------------------------------------------------------------*/
  /* ⑦複素数の絶対値自乗を取る（DSP関数）*/
  /*------------------------------------------------------------------------*/
	SquareMagnitudeCplx( FFT_POINT >>1 , frc_FftData, fr_PowerSpec ); 
	i4g_Time[7] = TMR1;
  /*------------------------------------------------------------------------*/
  /* ⑧正規化する*/
  /*------------------------------------------------------------------------*/
	for (i =0 ; i < FFT_POINT >> 1; i++)
	{
		fr_PowerSpec1[i] = fr_PowerSpec[i] << 5;
	}
	i4g_Time[8] = TMR1;

}

結果

解析する入力データを以下に示します。
・1kHzのサイン波を16kHzでサンプリングした入力データ(青線)
・窓関数データ(赤線)
・窓関数を適用した後の入力データ(緑線)

FFT演算結果

演算結果を下記に示します。1kHzにピークが現れている事が確認できます。

実行時間

以下に、FFTの実行時間の計測結果を示します。初期化を除くFFT実行関数（Calc_FFT）の処理時間は、コンパイラオプション1で60.31μsecとなりました。
dsPIC33Cで同様(コンパイラオプション1)に計測した場合、約120μsecだったため、理論通り約2倍の速度でFFTが実行できていることが確認できました。

コンパイラオプション0
(μsec)		コンパイラオプション1　
(μsec)
①回転因子の初期化22.5722.56
②窓関数の初期化23.9623.93
③窓関数の適用2.512.51
④入力データを複素数の実数に代入16.743.9
⑤FFTの実行45.9345.92
⑥ビット逆順ソート5.554.74
⑦複素数の絶対値自乗を取る1.281.27
⑧正規化する6.51.97
③～⑧合計78.5160.31

注意：メモリアロケーションについて

FFTを実行する際、回転因子とスペクトルデータの変数が正しくメモリに配置されていないと、期待通りの結果が得られないことがあります。結論として、スペクトルデータはXRAMに、回転因子はYRAMに割り当てることで正常な結果が得られます。

定義1 (無指定)・・・NG

fractcomplex twiddleFactors[FFT_POINT];
fractcomplex frc_FftData[ FFT_POINT ];

定義2 (XRAMにスペクトルデータ、YRAMに回転因子) ,near指定・・・OK

fractcomplex twiddleFactors[FFT_POINT] __attribute__ ( (space(ymemory),near));
fractcomplex frc_FftData[ FFT_POINT ] __attribute__ ( (space(xmemory),near));

定義3 （XRAMにスペクトルデータ、YRAMに回転因子) ,far指定・・・NG

fractcomplex twiddleFactors[FFT_POINT] __attribute__ ( (space(ymemory),far));
fractcomplex frc_FftData[ FFT_POINT ] __attribute__ ( (space(xmemory),far));

定義4　（XRAMに回転因子、YRAMにスペクトルデータ)・・・NG

fractcomplex twiddleFactors[FFT_POINT] __attribute__ ((section (".xbss, bss, xmemory"), aligned (FFT_POINT)));
fractcomplex frc_FftData[ FFT_POINT ] __attribute__ ((section (".ybss, bss, ymemory"), aligned (FFT_POINT)));

定義5 （XRAMにスペクトルデータ、YRAMに回転因子) ・・・OK

fractcomplex twiddleFactors[FFT_POINT] __attribute__ ((section (".ybss, bss, ymemory"), aligned (FFT_POINT)));
fractcomplex frc_FftData[ FFT_POINT ] __attribute__ ((section (".xbss, bss, xmemory"), aligned (FFT_POINT)));

記事についての注意点

本記事は慎重に内容を検討し正確さに努めておりますが、内容に誤りがあったとしても、この記事を参考にして生じた損害等については一切の責任を負いません。

コメント

タイトルとURLをコピーしました