dsPIC33AKシリーズ⑲ DSPの使い方4（LMSアルゴリズム）

dsPIC33Aシリーズ
2024.09.01

概要

今回は、dsPIC33AKシリーズのDSPを活用して、LMSアルゴリズムの検証を行います。

LMSアルゴリズムとは

LMS（Least Mean Square）アルゴリズムとは、最小2乗誤差を利用した適応フィルターです。目的の信号と実信号の差が最小になるように、係数を適応させるように動作します。アナログ回路では（実現する事は可能だが）現実的でない処理です。

なお今回の記事も以下の書籍（絶版）を参考にしております。

関連記事

記事リンク
第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)
第18回DSPの使い方3
(自己相関アルゴリズム)		dsPIC33AKシリーズ⑱ DSPの使い方3（自己相関アルゴリズム） – ぴくおの電子工作的な何かWP (electricpico.com)
第19回(本記事)DSPの使い方4
(LMSアルゴリズム)
(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)動作環境

動作説明

LMSアルゴリズムのブロック図

LMSアルゴリズムを用いたフィルタのブロック図を以下に示します。
FIRLMS関数は入力と出力の差分の2乗が最小になるようにをLMSアルゴリズムによって動的に導出され、FIRフィルタの係数を更新します。

システムブロック図

ソースコードのブロック図を下記に示します。

ソースコード

インクルードファイル

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

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

config.hダウンロード

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

Clock_Driver.cダウンロード

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

clock_driver.hダウンロード

ライブラリのインポート

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

(fig.3-2)ライブラリのインポート

CCP割り込みの中では以下のような処理を行っています。

①乱数（ノイズ）生成

信号周波数にノイズを加えるため、LFSRの値を使用しています。LFSR（Linear Feedback Shift Register) は読み出し毎に15ビットの擬似ランダム値を生成します。もともとスイッチング電源向けなどにノイズの高調波を拡散する目的などで使用されます。15ビットのLFSR値全体は自己相関特性を持ち、擬似ランダム・ノイズの発生源として望ましくないアプリケーションもあります。

(fig.3-3)LFSRブロック (データシート 70005539B – ｐ795)より抜粋
volatile int32_t RNDA = (LFSR << 14);
volatile int32_t RNDB = (LFSR << 14);

②信号遅延

LMSアルゴリズムの遅延量は1サンプル以上にする必要が有ります。この遅延量は環境によって左右されますので、適宜決めるようにします。なおVectorCopy関数でディレイラインをプッシュさせることも可能ですが、数サンプルであれば直接遅延処理を行った方が処理が速いです。

fl_InDelay[5] = fl_InDelay[4];
fl_InDelay[4] = fl_InDelay[3];
fl_InDelay[3] = fl_InDelay[2];
fl_InDelay[2] = fl_InDelay[1];
fl_InDelay[1] = fl_InDelay[0];
fl_InDelay[0] = sinData[(uint16_t)sinPointer] + 
             (int32_t)(LFSR << 14) - (int32_t)(LFSR << 14);
srcSamps = fl_InDelay[5];
refSamps = fl_InDelay[0];

③次回参照アドレスの算出と信号周波数の可変

SinRomから正弦波の瞬時値を取得する際には、sinPointerの値に基づいて参照アドレスを切り替えます。このsinPointerは、addPointerの値が加算されて更新されるため、addPointerの値を調整することで、出力される信号の周波数を制御することができます。言い換えれば、addPointerを増やせば周波数が高くなり、減らせば周波数が低くなる仕組みです。

sinPointer += addPointer;
if (sinPointer >= 1024)
{
    sinPointer = 0;
    Freq_side_Cnt++;
    if (Freq_side_Cnt > 5)
    {
	 Freq_side_Cnt = 0;
    }
}

④適応フィルタの実行

FIRLMS関数で適応フィルタの実行を行います。
係数は小さすぎると収束に時間がかかり、大きすぎると発散しやすくなります。環境によって左右されますので、適宜決めるようにします。

引数順序項目値
1実行回数1
2出力ポインタ&dstSamps
3遅延させた最後の入力データ&srcSamps
4FIRフィルタ構造体ポインタ&firFilter
5最新の入力データ& refSamps
6係数Q31(0.1)
FIRLMS ( 1,&dstSamps , &srcSamps,&firFilter , & refSamps ,Q31(0.1));

⑤DMA送信

入力信号と出力信号をPCに出力するため、sprintf関数でバッファに格納した後に、DMA転送を実行します

DMA0CNT = sprintf(TxBuf,"%d,%d\r\n",refSamps,dstSamps) ;
DMA0SRC = (unsigned int)& TxBuf; //転送元
DMA0CHbits.CHEN = 1U;  //DMA0チャネル有効化
DMA0CHbits.CHREQ = 1U; //DMA0ソフトウェア転送開始

ソースコード全体

#include <stdio.h>
#include <stdlib.h>
#include <xc.h>
#include "config.h"
#include "clock_driver.h"
#include <dsp.h>
/*----------------------------------------------------------------------------*/
/* 定数定義*/
/*----------------------------------------------------------------------------*/
#define TXBUF_SIZE 256
#define FIR_TAP 100
#define FSAMP   1000	//サンプリング周波数
#define FSIGA    12	//信号周波数A
#define FSIGB    25	//信号周波数B
/*----------------------------------------------------------------------------*/
/* グローバル変数定義*/
/*----------------------------------------------------------------------------*/
fractional fl_InDelay[8];
fractional srcSamps;
fractional dstSamps;
fractional refSamps;
fractional sinData[1024];
float sinPointer;
float addPointer;
int Freq_side_Cnt;
int8_t TxBuf[TXBUF_SIZE];

fractional firDelay[FIR_TAP] __attribute__ ((space(ymemory), near));
fractional firCoeff[FIR_TAP] __attribute__ ((space(xmemory), near));
FIRStruct firFilter;

/*----------------------------------------------------------------------------*/
/**
* @fn          main(int argc, char** argv) 
* @brief       PROJECT_4_1_4_DSP_FIRLMS
* @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) 
{
    /*------------------------------------------------------------------------*/
    /*　クロック初期化*/
    /*------------------------------------------------------------------------*/
        vdg_Clock_Set_Register();    /* クロック初期化 */
  　/*-----------------------------------------------------------------------*/
    /* ピン設定*/
    /*-----------------------------------------------------------------------*/
	RPOR12bits.RP50R = 9u;		//UART1出力 RP50
	RPINR9bits.U1RXR = 52u;		//UART1入力 RP52
    /*-----------------------------------------------------------------------*/
    /* UART初期化*/
    /*-----------------------------------------------------------------------*/
	U1CON = 0x00000000u;
	U1STAT = 0x00000000u;
	U1BRG = (uint32_t) (100000000 / 460800);
	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;

	INTCON4bits.OVATE = 0;
		
	for (int u2l_Loop = 0;u2l_Loop < 1024 ;u2l_Loop ++)
        {
	  sinData[u2l_Loop ] = Q31(sin((float)u2l_Loop /1024 * 2 * PI)) >> 2;
        }
    /*-----------------------------------------------------------------------*/
    /* DMA設定 */
    /*-----------------------------------------------------------------------*/
	DMACONbits.ON = 1U;        //DMA有効
	DMACONbits.PRIORITY = 0U;  //Fixed
	DMACONbits.SIDL = 0UL;
	DMABUF = 0x00000000u;
	DMALOW = 0x4000U;
	DMAHIGH = 0x7FFFU;

	DMA0CH = 0x00000000U;
	DMA0CHbits.SIZE = 0U;	//8Bit
	DMA0CHbits.TRMODE = 0U; //ワンショット
	DMA0CHbits.DAMODE = 0U; //
	DMA0CHbits.SAMODE = 1U; //Increment

	DMA0SEL = 0x10U;
	DMA0STAT = 0x00000000U;
	DMA0SRC = (unsigned int)& TxBuf; //転送元
	DMA0DST = (unsigned int)& U1TXB; //転送先
	DMA0CNT = 8u;                           //転送回数
	DMA0CLR = 0x00000000U;
	DMA0SET = 0x00000000U;
	DMA0INV = 0x00000000U;
	DMA0MSK = 0x00000000U;
	DMA0PAT = 0x00000000U;

    /*------------------------------------------------------------------------*/
    /* Filter初期化 */
    /*------------------------------------------------------------------------*/
        FIRStructInit( &firFilter,          // フィルター構造体
                        FIR_TAP,            // フィルター係数タップ数
                        firCoeff,	    // フィルタ係数ポインタ
                        0,		    // フィルタ係数PSVページ
                        firDelay);          // ディレイバッファ
	FIRDelayInit( &firFilter );
		
	for(int i = 0; i < 100; i++) {firCoeff[i] = 0;}
　　　　addPointer = 1024/(FSAMP / FSIGA);
    /*-----------------------------------------------------------------------*/
    /* SCCP1設定 */
    /*-----------------------------------------------------------------------*/
	CCP1CON1 = 0x00000000u;
	CCP1CON1bits.T32 = 1u;
	CCP1CON2 = 0x00000000u;
	CCP1CON3 = 0x00000000u;
	CCP1PR = 100000000 / FSAMP;
		
	IPC6bits.CCT1IP = 5;
	IFS1bits.CCT1IF = 0;
	IEC1bits.CCT1IE = 1u;
	CCP1CON1bits.ON = 1u;
		
	INTCON1bits.GIE = 1;	
		
    /*------------------------------------------------------------------------*/
    /* メインルーチン */
    /*------------------------------------------------------------------------*/
	while (1)
	{
	    if (Freq_side_Cnt < 3)
	    {
		addPointer = 1024/(FSAMP / FSIGA);
	    }
	    else
	    {
		addPointer = 1024/(FSAMP / FSIGB);
	    }
	}
        return EXIT_SUCCESS;
}


void __attribute__((interrupt, no_auto_psv)) _CCT1Interrupt(void)
{
    /*------------------------------------------------------------------------*/    
    /* ①乱数生成 */
    /*------------------------------------------------------------------------*/
	volatile int32_t RNDA = (LFSR << 14);
	volatile int32_t RNDB = (LFSR << 14);
    /*------------------------------------------------------------------------*/
    /* ②信号の遅延 */
    /*------------------------------------------------------------------------*/
	fl_InDelay[5] = fl_InDelay[4];
	fl_InDelay[4] = fl_InDelay[3];
	fl_InDelay[3] = fl_InDelay[2];
	fl_InDelay[2] = fl_InDelay[1];
	fl_InDelay[1] = fl_InDelay[0];
	fl_InDelay[0] = sinData[(uint16_t)sinPointer] + 
             (int32_t)(LFSR << 14) - (int32_t)(LFSR << 14);
	srcSamps = fl_InDelay[5];
	refSamps = fl_InDelay[0];
    /*------------------------------------------------------------------------*/
    /* ③次回参照アドレスの算出と信号周波数の可変 */
    /*------------------------------------------------------------------------*/
	sinPointer += addPointer;
	if (sinPointer >= 1024)
	{
	    sinPointer = 0;
	    Freq_side_Cnt++;
	    if (Freq_side_Cnt > 5)
	    {
		 Freq_side_Cnt = 0;
	    }
	}
    /*------------------------------------------------------------------------*/
    /* ④適応フィルタの実行 */
    /*------------------------------------------------------------------------*/
	FIRLMS ( 1,&dstSamps , &srcSamps,&firFilter , & refSamps ,Q31(0.1));
    /*------------------------------------------------------------------------*/
    /* ⑤DMA送信 */
    /*------------------------------------------------------------------------*/
        DMA0CNT = sprintf(TxBuf,"%d,%d\r\n",refSamps,dstSamps) ;
	DMA0SRC = (unsigned int)& TxBuf; //転送元
	DMA0CHbits.CHEN = 1U;  //DMA0チャネル有効化
	DMA0CHbits.CHREQ = 1U; //DMA0ソフトウェア転送開始

	IFS1bits.CCT1IF = 0;
}

結果1

サンプリング周波数等を以下の値に設定した時の入力信号と出力信号を下記に示します。

項目設定値
サンプリング周波数(Hz)1000
遅延量5
信号周波数A(Hz)12
(fig.3-4)設定値

■係数 = Q31(0.1)

■係数 = Q31(0.25)

■係数 = Q31(0.5)

結果2

サンプリング周波数等を以下の値に設定した時の入力信号と出力信号を下記に示します。

項目設定値
サンプリング周波数(Hz)1000
遅延量5
信号周波数A(Hz)12
信号周波数B(Hz)25
係数Q(0.1)
(fig.3-8)設定値

記事についての注意点

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

コメント

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