概要
dsPICシリーズは「dsP」と頭文字がつくように、DSP(*1)が内蔵されていることは前回の記事で紹介しました。今回はこのDSPを用いFFT(*2)を構成します。
なお本記事は絶版書になりますがCQ出版社から出版されていた以下の著書を参考にしております。
[絶版2017.1.24] dsPIC基板で始めるディジタル信号処理 (cqpub.co.jp)
FFTとは高速フーリエ変換(Fast Fourier Transform)の略。フーリエ変換とは時間領域の関数を周波数領域の関数に変換することを指し、それを高速に演算するアルゴリズムの事
ここでは詳しいアルゴリズムを紹介できませんしませんが、どのようなものなのかを紹介するサイト様は調べればたくさん出てきますので、いくつかご紹介します。
FFT（高速フーリエ変換）を完全に理解する話
FFT解析とは？【原理とグラフの見方を3ステップでわかりやすく解説】
【フーリエ解析05】高速フーリエ変換(FFT)とは？内側のアルゴリズムを解説！【解説動画付き】
3．周波数分析の基礎 – Oita University
いかがでしたでしょうか？
ちなみに私は「時系列データにどの周波数成分がどれ位含まれるかを高速に算出するアルゴリズム」位しか理解していません。
ハードウェア構成と制御ブロック
今回のプログラムはハードウェアは不要でシミュレータ上で動作させます。
レジスタ
レジスタは未使用です。
ソースコード
DSP命令を利用する時は(fig.12)のようにプロジェクトのLibrariesに「libdsp-elf.a」をインポートします。
また「dsp.h」をインクルードします。
#include <dsp.h>DSPによるFFTの実行ですが、今回は予め ダミーデータとして1kHzのサイン波を16kHzでサンプリングした128ポイントのデータをfr_InputDataに格納します。
//入力データ 128サンプル 1000Hz/16kHz fs
fractional fr_InputData[ FFT_POINT ] =
{
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540 ,
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540 ,
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540 ,
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540 ,
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540 ,
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540 ,
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540 ,
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540
};初期設定として以下のコードのように、回転因子の初期化と窓関数の初期化を行います。
//回転因子の初期化
tFactor = TwidFactorInit(FFT_POWEROF2,(fractcomplex *)twiddleFactors,0);
//窓関数の初期化
HanningInit(FFT_POINT,fr_Window);初期化が終わった時点でFFTを実行します
//FFTの実行
Calc_FFT(&fr_InputData,FFT_POINT,FFT_POWEROF2);コンフィグレーション設定についてはコンフィグレーション設定に記載しております。
コピーして下記のソースコードの「 //ここにコンフィグレーション設定を挿入する// 」の位置に挿入してください。
クロック設定用関数 vds_Main_Init_Clock_Register(); のソースコードはクロック設定のページに記載しております。
コピーして下記のソースコードの「 //ここにクロック設定ソースをコピペする// 」の位置に挿入してください。
/*----------------------------------------------------------------------------*/
/* <Chapter> CHAPTER_4_1_3_DSP_FFT */
/* <Function> DSPを利用したFFT */
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
/* コンフィグレーション設定*/
/*----------------------------------------------------------------------------*/
//ここにコンフィグレーション設定を挿入する//
/*----------------------------------------------------------------------------*/
/* インクルードファイル*/
/*----------------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <xc.h>
#include <libpic30.h>
#include <dsp.h>
/*----------------------------------------------------------------------------*/
/* 定数定義*/
/*----------------------------------------------------------------------------*/
#define FFT_POWEROF2 7u
#define FFT_POINT ( 1 << FFT_POWEROF2)
/*----------------------------------------------------------------------------*/
/* 変数定義*/
/*----------------------------------------------------------------------------*/
//回転因子
fractcomplex twiddleFactors[FFT_POINT] __attribute__ ((section (".xbss, bss, xmemory"), aligned (FFT_POINT)));
//回転因子ポインタ
fractcomplex *tFactor = (fractcomplex * )&twiddleFactors;
//入力データ 128サンプル 1000Hz/16kHz fs
fractional fr_InputData[ FFT_POINT ] =
{
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540 ,
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540 ,
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540 ,
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540 ,
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540 ,
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540 ,
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540 ,
0 ,12539 ,23169 ,30272 ,32767 ,30272 ,23169 ,12539 ,
0 ,-12540 ,-23170 ,-30273 ,-32767 ,-30273 ,-23170 ,-12540
};
// スペクトル(複素数)
fractcomplex frc_FftData[ FFT_POINT ] __attribute__ ( (space(ymemory),far));
// パワースペクトル（実数）
fractional fr_PowerSpec[ FFT_POINT/2 ];
fractional fr_PowerSpec1[ FFT_POINT/2 ];
// 窓関数用
fractional fr_Window[ FFT_POINT ];
/*----------------------------------------------------------------------------*/
/* クロック設定 */
/*----------------------------------------------------------------------------*/
//ここにクロック設定ソースをコピペする//
/*----------------------------------------------------------------------------*/
/**
* @fn main(int argc, char** argv)
* @brief CHAPTER_4_1_1_DSP_FIR
* @param[in] argc argument count
* @param[in] argv argument vector
* @retval EXIT_SUCCESS 成功
* @retval EXIT_FAILURE 失敗
* @details DSPによるFFT
* @note
*/
/*----------------------------------------------------------------------------*/
int main(int argc, char** argv)
{
//回転因子の初期化
tFactor = TwidFactorInit(
FFT_POWEROF2,
(fractcomplex *)twiddleFactors,
0);
//窓関数の初期化
HanningInit(FFT_POINT,fr_Window);
//FFTの実行
Calc_FFT(&fr_InputData,FFT_POINT,FFT_POWEROF2);
while (1)
{
}
return EXIT_SUCCESS;
}
/*----------------------------------------------------------------------------*/
/**
* @fn Calc_FFT(_s2 *s2a_Data ,_u2 u2a_FFT_Points,_u2 u2a_FFT_PO2)
* @brief FFTの実行
* @param[in] *s2a_Data 入力データ配列
* @param[in] u2a_FFT_Points FFTポイント数
* @param[in] u2a_FFT_PO2 FFTポイント数 2＾ｎ
* @details FFT実行
* @note
*/
/*----------------------------------------------------------------------------*/
void Calc_FFT(_s2 *s2a_Data ,_u2 u2a_FFT_Points,_u2 u2a_FFT_PO2)
{
int i;
//窓関数の適用（DSP関数）[548サイクル]
VectorWindow(FFT_POINT,
(fractional *)s2a_Data,
(fractional *)s2a_Data,
(fractional *)fr_Window);
//入力データを複素数の実数に代入 [1283サイクル]
for (i = 0; i < u2a_FFT_Points;i++)
{
frc_FftData[i].real = *s2a_Data++;
frc_FftData[i].imag = 0;
}
// FFTの実行（DSP関数）[8765サイクル]
FFTComplexIP(
u2a_FFT_PO2,
(fractcomplex *)frc_FftData,
(fractcomplex *)twiddleFactors,
COEFFS_IN_DATA);
// ビット逆順ソート（DSP関数）[1097サイクル]
BitReverseComplex( u2a_FFT_PO2, frc_FftData );
// 複素数の絶対値自乗を取る（DSP関数）[290サイクル]
SquareMagnitudeCplx( FFT_POINT/2, frc_FftData, fr_PowerSpec );
}結果
解析する入力データを(fig.2)に示します。
・1kHzのサイン波を16kHzでサンプリングした入力データ(fir_InputData/青線)
・窓関数データ(fir_WindowData/赤線)
・窓関数を適用した後の入力データ(fir_InputData/緑線)
窓関数を適用することで両端のデータをスムーズにつなげます。
窓関数を用いる詳細は以下のサイト様を参考下さい。
窓関数を用いる理由
FFT実行後のスペクトルは周波数順に並んでいないため、ビットリバース処理を行います。
ビットリバース処理後の複素数データを(fig.3)に示します。
複素数の実部と虚部の2乗を加算した値の平方根が振幅スペクトルAになります。
パワーは振幅の2乗に比例します。つまり平方根が取れますので実部と虚部の2乗を加算した値がパワースペクトルとなります。
DSP命令のSquareMagnitudeCplx関数は複素数の実部2乗と虚部の2乗を加算して出力しますので、直接パワースペクトルが求められます。
SquareMagnitudeCplx( FFT_POINT/2, frc_FftData, fr_PowerSpec ); SquareMagnitudeCplx関数で演算したパワースペクトル(fr_PowerSpec)を2倍で正規化したデータを(fig.5)に示します。1kHzにピークが現れている事が確認できます。
また窓関数を補正した値を補正した値を(fig.6)に示します。
なお初期化以外のFFT実行関数（Calc_FFT)の処理サイクルは12,004サイクルとなりました。（XC16 2.00/コンパイラオプション1での結果)
100MHzで実行したとすると所要サイクルは約120usecとなります。
Microchipの16bit言語ツールというガイドにDSPライブラリの説明や実行サイクル数の記載が有りますが、どのDSPライブラリも実行サイクルは理想サイクルよりも数％多くかかっています。理由は不明ですが、これより早く実行させる場合はC言語ライブラリでは無くアセンブラを使うしか無いと思います。
| DSP関数 | 説明 | 理想サイクル数 | 実行サイクル |
| VectorMultiply | VectorMultiply は、ソース1 ベクタ内の各要素値とソース2 ベク タ内の対応する各要素値の積をとり、結果をディステネーション・ベ クタの対応する各要素に格納します。 | 17＋ 4(NUM) = 529 | – |
| VectorWindow | VectorWindow は、ウインドウを指定されたソース・ベクタに適用し、ウインドウ化された結果のベクタをディステネーション・ベクタ へ格納します。 | 3 + VectorVectorMultiply =532 | 548 |
| FFTComplexIP | FFTComplexIP は、ソース複素数ベクタの離散的フーリェ変換をイン・プレイス計算します。 | 8,485 調整係数がX- メモリの場合 | 8,765 |
| BitReverseComplex | BitReverseComplex は、複素数ベクタの要素のビット順を逆にします。 | 945 | 1,097 |
| SquareMagnitudeCplx | SquareMagnitudeCplx は、複素数ソース・ベクタ内の各要素の大きさの2 乗を計算します。 | 20 + 3(numElems) = 212 | 290 |
コメント