概要
今回は、dsPIC33AKシリーズのDSPを活用して、自己相関アルゴリズムの検証を行います。
自己相関（Auto-correlation）とは、時系列データにおいて、データの過去の値と現在の値との間にどのような関係があるかを分析する手法です。自己相関を調べることで、データに周期性やトレンドがあるか、ノイズなのかを判断するのに役立ちます
なお今回の記事も以下の書籍（絶版）を参考にしております。
関連記事
| 関連記事 | リンク |
| DSPライブラリの使い方1（デジタルフィルタ） | DSPライブラリの使い方1（デジタルフィルタ） – ぴくおの電子工作的な何かWP (electricpico.com) |
| DSPライブラリの使い方2（FFT） | DSPライブラリの使い方2（FFT） – ぴくおの電子工作的な何かWP (electricpico.com) |
開発環境
開発環境を以下に示します。
今回はCPUボード(EV02G02A)で動作確認を行います。
| 項目 | 値 | リンク |
| ベースボード | dsPIC33A CURIOSITY PLATFORM DEVELOPMENT BOARD | dsPIC33A Curiosity Platform Development Board User’s Guide (microchip.com) |
| CPUボード(EV68M17A) | EV68M17A – dsPIC33AK128MC106 Motor Control DIM | dsPIC33AK128MC106 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.20 | MPLAB® X IDE | Microchip Technology |
| コンパイラ | MPLAB XC DSC v3.10 | MPLAB® XC DSC Compiler | Microchip Technology |
動作説明
自己相関は、その名の通り、ある時系列データが時間的にどの程度自身と相関しているか（つまり、データが周期的であるかどうか）を測定するために使用されます。これは、以下の式で表されます。
ACR[n] :ラグ nnn における自己相関関数の値を表します。ラグ nnn は、データがどれだけ前の時間点にシフトしているかを表します。
x[k]: 時系列データの k 番目の値を表します。
x[k+n]: 時系列データの k+n 番目の値を表します。
つまりn が変化することで、データを時間的にシフトさせた場合にその自己相関を計算します。
自己相関アルゴリズムのメリット・デメリット
メリット:
- 精度: ピークの位置を特定することで、信号の基本周波数を非常に高い精度で求めることができます。
- 周期的信号の特定: 特に周期的な信号の基本周波数を求めるのに有効です。音声処理や振動解析でよく使われます。
- 信号長に依存しない: FFTのように、信号の長さに依存せず、短い信号でも有効です。
デメリット:
- 計算量: 自己相関はO(N^2)の計算量を持つため、長いデータセットに対しては計算が遅くなります。
- ノイズの影響: ノイズの影響を受けやすく、信号にノイズが含まれると精度が低下することがあります。
- 多重周波数: 基本周波数以外の複数の周波数成分が含まれている場合、解析が困難になることがあります。
ソースコード
インクルードファイル
コンフィグレーションファイル、クロック設定ファイルは以下のファイルをインクルードしてください。
■コンフィグレーションファイル (config.h)
■クロック設定ソースファイル(Clock_Driver.c)
■クロックヘッダーファイル(clock_driver.h)
ライブラリのインポート
DSPを使用する場合、dsp.hのインクルードとlibdsp-elf.aのインポートが必須です。
VectorDotProduct関数について
今回使用するVectorDotProduct関数は積和演算を行ってくれる便利な関数です。
| 引数順 | 値 | |
| 1 | 積和回数 | 64 |
| 2 | 第一変数先頭ポインタ | fl_Data |
| 3 | 第二変数先頭ポインタ | fl_Data + i |
fl_ACR[i] = VectorDotProduct(64,fl_Data,fl_Data + i);ただし積和演算の特性上、アキュムレータが飽和/オーバーフローする可能性が有ります。これについては後述します。
ソースコード全体
#include <stdio.h>
#include <stdlib.h>
#include <xc.h>
#include "config.h"
#include "clock_driver.h"
#include <dsp.h>
/*----------------------------------------------------------------------------*/
/* 定数定義*/
/*----------------------------------------------------------------------------*/
#define FREQ_A 8000.0
#define FREQ_B 000.0
#define FREQ_C 12000.0
#define GAIN_A 0.9
#define GAIN_B 0.5
#define GAIN_C 0.0
#define SFREQ 48000
/*----------------------------------------------------------------------------*/
/*変数定義*/
/*----------------------------------------------------------------------------*/
fractional fl_ACR[128];
fractional fl_Data[128];
static void Uart_Send_Blocking_TX_Data(uint8_t *u1pa_Data);
/*----------------------------------------------------------------------------*/
/**
* @fn main(int argc, char** argv)
* @brief PROJECT_1_3_3_DSP_ACR
* @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];
float f4_Time,f4_DummyA,f4_DummyB,f4_DummyC;
/*------------------------------------------------------------------------*/
/* クロック初期化*/
/*------------------------------------------------------------------------*/
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;
for (int u2l_Loop = 0;u2l_Loop < 128 ;u2l_Loop ++)
{
f4_Time = (float)u2l_Loop/(SFREQ);
f4_DummyA = sin(f4_Time * FREQ_A * 2 * PI) * GAIN_A ;
f4_DummyB = sin(f4_Time * FREQ_B * 2 * PI) * GAIN_B ;
f4_DummyC = sin(f4_Time * FREQ_C * 2 * PI) * GAIN_C ;
fl_Data[u2l_Loop] = Q31(f4_DummyA + f4_DummyB + f4_DummyC) >> 4;
}
/*------------------------------------------------------------------------*/
/* 自己相関を求める */
/*------------------------------------------------------------------------*/
for (int i = 0;i< 32;i ++)
{
fl_ACR[i] = VectorDotProduct(64,fl_Data,fl_Data + i);
}
/*------------------------------------------------------------------------*/
/* 結果の出力 */
/*------------------------------------------------------------------------*/
for (int i = 0;i< 128;i ++)
{
sprintf(TxTemp,"%d,%d,%d\r\n",i,fl_Data[i],fl_ACR[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 ++;
}
}
結果1
周波数とゲインを以下の値に設定した時の入力信号と出力(相関)信号を下記に示します。
| A | B | C | |
| FREQ | 2000 | 4000 | 8000 |
| GAIN | 0.0625 | 0 | 0 |
出力信号を要素順にプロットすると、24番目のデータが極大となっていることが確認できました。これは、2kHzの信号を48kHzでサンプリングしているため、48/2 = 24 サンプルで位相が1回転します。このように、自己相関アルゴリズムでは位相がそろったところで極大値を持ちます。
結果2
周波数とゲインを以下の値にした時の入力信号と出力(相関)信号を下記に示します。
| A | B | C | |
| FREQ | 2000 | 5000 | 8000 |
| GAIN | 0.0625 | 0.03125 | 0 |
結果1のように簡単に算出できる一方、デメリットの項目にも記載したとおり、基本周波数以外の周波数が重畳している場合、解析が困難になることがあります。
結果3
周波数とゲインを以下の値にした時の入力信号と出力(相関)信号を下記に示します。
| A | B | C | |
| FREQ | 2000 | 4000 | 8000 |
| GAIN | 0.25 | 0 | 0 |
出力データが-1 ～ 1 で飽和してしまいました。これは積和演算の特性上、避けられない現象です。
アキュムレータ自体は72bit幅ですので飽和はまだ発生していません（SRbits.OA = 0)が、戻り値は丸めロジックで32bit幅になりますので、このライトバックの時点で飽和が発生します（SRbits.OAB = 1)
ちなみに、dsPIC33Aシリーズに関してまだ理解できていない部分も多々ありますので、誤解があるかもしれませんが、dsPIC33AシリーズではSRレジスタやワーキングレジスタがメモリマップされていないため、直接SRレジスタの値を取得して演算が飽和したかどうかを確認することはできないようです。INTCON4bits.OVATEフラグをセットすれば算術演算トラップが発生しますが、それ以外の方法で飽和を確認する手段については現在調査中です。
記事についての注意点
本記事は慎重に内容を検討し正確さに努めておりますが、内容に誤りがあったとしても、この記事を参考にして生じた損害等については一切の責任を負いません。
コメント