FFTComplexIP関数,IFFTComplexIP 5月号(PP226-229)を使います。
少々トラ技から逸脱しますが、信号として16QAMのOFDM信号(もどき)を扱います。
dsPICで処理するための信号を作ります。ここでは16QAMの信号としました。
16QAMデータ ⇒ IQデータマッピング ⇒IFFT/FFTポイント位置に並べる
⇒IFFTする ---- ⇒ FFTする ⇒IQデータに戻る ⇒16QAMデータに戻る
(1)16QAMの信号を作成
計算ソフトoctaveを使ってコンスタレーションを作成してみます
M=16;%16QAM
x = [0:M-1];
y2 = qammod(x,M);
%Matlabと違ってmode(gray codeやbinary)の指定ができませんのでbinaryのまま計算です
scatterplot(y2,1,0,'b.'); % スケールされた信号点配置をプロット
legend 16QAM;
title('Constellation for Binary-Coded 16-QAM');
% 点に番号を付けるテキスト注釈を含む
hold on; % 注釈を同じ図に置く
annot = dec2bin([0:length(y2)-1],log2(M));
text(real(y2)+0.15,imag(y2),annot);
hold off;
計算結果はx=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]に対して
y2=[-3+3i,-3+1i,-3-1i,-3-3i,-1+3i,-1+1i,-1-1i,-1-3i,1+3i,1+1i,1-1i,1-3i,3+3i,3+1i,3-1i,3-3i]
%FFT数を64ポイントにして、オリジナルデータのスペクトラムを表示ですFig.2
yy2=[y2 y2 y2 y2]; %64 points(データへの0挿入はしてません)
fpoints = 16*2*2; %64points
figure(1);
n=1:fpoints;
subplot(2,1,1)
yt(1,:) = abs(yy2);
yt(2,:) = real(yy2);
yt(3,:) = imag(yy2);
plot(n,yt) % for check
lgend ('abs','I','Q');
title('yy2')
%次にOctaveに戻って IFFTの計算です
I,Q成分も表示してみます。(Fig.3)
figure(2);
n=1:fpoints;
subplot(2,1,1)
yift= ifft(yy2);
yit(1,:) = abs(yift);
yit(2,:) = real(yift);
yit(3,:) = imag(yift);
plot(n,yit) % plot
legend ('abs','I','Q');
title('yy2-ifft')
さらに、fftして元に戻してみます。(Fig.4)
n=1:fpoints;
subplot(2,1,2)
yift_ft= fft(yift);
yift(1,:)= abs(yift_ft);
yift(2,:) = real(yift_ft);
yift(3,:) = imag(yift_ft);
plot(n,yift); % for checkf
legend ('abs','I','Q');
title('yy2-ifft-fft');
Fig.2とFig.4が同じくなりました。IFFTしてさらにFFTすることでデータが元に戻りました。
これと同じことを以下dsPICで行ってみます。
octaveと同じことを行ってみます。
CDに付属のサンプルソースCE018_FFT_DSPLibを触りながら進めます
関係しそうなファイルは
inputsignal_square1khz.c
main_FFTExample.c
twiddleFactors.c
fft.h
linkerscript.gld
inputsignal_square1khz.cに固定小数点に変換したIQデータをおきます。
元データ:y2=[-3+3i,-3+1i,-3-1i,-3-3i,-1+3i,-1+1i,-1-1i,-1-3i,1+3i,1+1i,1-1i,1-3i,3+3i,3+1i,3-1i,3-3i]
変換したデータ
fractcomplex sigCmpx[FFT_BLOCK_LENGTH] __attribute__ ((section (".ydata, data, ymemory"),
aligned (FFT_BLOCK_LENGTH * 2 *2))) =
{
0xC000,0x4000,0xC000,0x1555,0xC000,0xEAAB,0xC000,0xC000,
0xEAAB,0x4000,0xEAAB,0x1555,0xEAAB,0xEAAB,0xEAAB,0xC000,
0x1555,0x4000,0x1555,0x1555,0x1555,0xEAAB,0x1555,0xC000,
0x4000,0x4000,0x4000,0x1555,0x4000,0xEAAB,0x4000,0xC000,
0xC000,0x4000,0xC000,0x1555,0xC000,0xEAAB,0xC000,0xC000,
0xEAAB,0x4000,0xEAAB,0x1555,0xEAAB,0xEAAB,0xEAAB,0xC000,
0x1555,0x4000,0x1555,0x1555,0x1555,0xEAAB,0x1555,0xC000,
0x4000,0x4000,0x4000,0x1555,0x4000,0xEAAB,0x4000,0xC000,
0xC000,0x4000,0xC000,0x1555,0xC000,0xEAAB,0xC000,0xC000,
0xEAAB,0x4000,0xEAAB,0x1555,0xEAAB,0xEAAB,0xEAAB,0xC000,
0x1555,0x4000,0x1555,0x1555,0x1555,0xEAAB,0x1555,0xC000,
0x4000,0x4000,0x4000,0x1555,0x4000,0xEAAB,0x4000,0xC000,
0xC000,0x4000,0xC000,0x1555,0xC000,0xEAAB,0xC000,0xC000,
0xEAAB,0x4000,0xEAAB,0x1555,0xEAAB,0xEAAB,0xEAAB,0xC000,
0x1555,0x4000,0x1555,0x1555,0x1555,0xEAAB,0x1555,0xC000,
0x4000,0x4000,0x4000,0x1555,0x4000,0xEAAB,0x4000,0xC000
};
ここでfractcomplexは小数複素数ベクタで各ベクタは対であらわされReal(2バイト)、Imag(2バイト)の
順に格納されるます。このデータは1ビットの符号と15ビットの小数であらわされるため”1.15データ”と呼ばれ、
±1の範囲の値しか表現できませんが、FFT関数では[-0.5,+0.5]の範囲のデータを扱います。
また、IQデータを固定小数点に変換する際にはofficeのExcelが便利です。参考文献(4)に紹介されてます。
DSP librayには
struct {
fractional real;
fractional imag;
}fractcomplex;
が定義されています。
次に
main_FFTExample.cでIFFTを演算します
TwidFactorInit (LOG2_BLOCK_LENGTH, &twiddleFactIFFT2[0], 1);
IFFTComplexIP(LOG2_BLOCK_LENGTH,&sigCmpx[0], &twiddleFactIFFT2[0], COEFFS_IN_DATA);
さらにFFTして16QAMの信号に戻します。
FFTComplexIP (LOG2_BLOCK_LENGTH, &sigCmpx[0], (fractcomplex *) __builtin_psvoffset(&twiddleFactors[0]), (int) __builtin_psvpage(&twiddleFactors[0]));
この結果を表計算ソフトExelでコンスタレーションの比較を行います。
予定通り元に戻りました。正常に動作していることがわかります。
FFT、IFFTをつかうとき、Microchipのドキュメントにバグがあるようです。(参考文献5)
BitReverseComplex関数の位置に気をつけてください。
// *********************以下プログラムです*************************************
// ************************* main_FFTExample.c********************************
#include < p30F2012.h >
#include < dsp.h >
#include "fft.h"
#include < uart.h >
#include < stdio.h >
#include < string.h >
// UART control register
unsigned int _U1BRG;
unsigned int _U1MODE;
unsigned int _U1STA;
#define Fcy 7370000*4
#define BaudRate 115200
// data
unsigned int ResultData;
/* Device configuration register macros for building the hex file */
_FOSC(CSW_FSCM_OFF & XT_PLL8); /* XT with 8xPLL oscillator, Failsafe clock off */
_FWDT(WDT_OFF); /* Watchdog timer disabled */
_FBORPOR(PBOR_OFF & MCLR_EN); /* Brown-out reset disabled, MCLR reset enabled */
_FGS(CODE_PROT_OFF); /* Code protect disabled */
fractcomplex twiddleFactIFFT2[FFT_BLOCK_LENGTH/2] __attribute__ ((space(xmemory)));
extern fractcomplex sigCmpx[FFT_BLOCK_LENGTH] /* Typically, the input signal to an FFT */
__attribute__ ((section (".ydata, data, ymemory"), /* routine is a complex array containing samples */
aligned (FFT_BLOCK_LENGTH * 2 *2))); /* of an input signal. For this example, */
extern const fractcomplex twiddleFactors[FFT_BLOCK_LENGTH/2] /* Twiddle Factor array in Program memory */
__attribute__ ((space(auto_psv), aligned (FFT_BLOCK_LENGTH*2)));
int main(void)
{
int i = 0;
int n;
char buffer_in[20];
TRISCbits.TRISC13=0; // RC13 is outport for led
PORTCbits.RC13=0; // green LED on (RC13=1)
// Use Alternate Interrupt Vector Table
INTCON2bits.ALTIVT=1;
// UART1モジュールのオフを確認します
CloseUART1();
// UART1割り込みを設定します
ConfigIntUART1( UART_RX_INT_DIS & UART_TX_INT_DIS);
// UART1を初期化します
// BaudRate=Fcy/(16(U1BRG+1))
// =29.48MHz/(16(15+1))
// =115156
// Error=(115156-115200)/115200
// =-0.0382%
_U1BRG= ((Fcy/BaudRate)/16);
// alternate I/O (U1ATX/U1ARX), N81
_U1MODE= UART_EN &
UART_IDLE_CON &
UART_ALTRX_ALTTX &
UART_DIS_WAKE &
UART_DIS_LOOPBACK &
UART_EN_ABAUD &
UART_NO_PAR_8BIT &
UART_1STOPBIT;
_U1STA= UART_INT_TX_BUF_EMPTY &
UART_TX_PIN_NORMAL &
UART_TX_ENABLE &
UART_INT_RX_CHAR &
UART_ADR_DETECT_DIS &
UART_RX_OVERRUN_CLEAR;
// Open UART1
OpenUART1(_U1MODE, _U1STA, _U1BRG);
strcpy(buffer_in,"Start\0");
putsUART1((unsigned int*)buffer_in);//送信する文字列をUART送信バッファに書き込む関数。NULL文字で終了する。
strcpy(buffer_in,"init,\0");
putsUART1((unsigned int*)buffer_in);//送信する文字列をUART送信バッファに書き込む関数。NULL文字で終了する。
/* Generate TwiddleFactor Coefficients */
TwidFactorInit (LOG2_BLOCK_LENGTH, &twiddleFactIFFT2[0], 1); /* We need to do this only once at start-up */
strcpy(buffer_in,"--twiddle--,\0");
putsUART1((unsigned int*)buffer_in);//送信する文字列をUART送信バッファに書き込む関数。NULL文字で終了する。
for ( n=0 ; n < FFT_BLOCK_LENGTH/2 ; n++ )
{
sprintf(buffer_in, "%X,\0", twiddleFactIFFT2[n].real);
putsUART1((unsigned int*)buffer_in);//送信する文字列をUART送信バッファに書き込む関数。NULL文字で終了する。
sprintf(buffer_in, "%X,\0\n", twiddleFactIFFT2[n].imag);
putsUART1((unsigned int*)buffer_in);//送信する文字列をUART送信バッファに書き込む関数。NULL文字で終了する。
}
strcpy(buffer_in,"--ifft--,\0");
putsUART1((unsigned int*)buffer_in);//送信する文字列をUART送信バッファに書き込む関数。NULL文字で終了する。
/* Perform IFFT operation */
IFFTComplexIP(LOG2_BLOCK_LENGTH,&sigCmpx[0], &twiddleFactIFFT2[0], COEFFS_IN_DATA);
for (n=0 ; n < FFT_BLOCK_LENGTH ;n++)
{
sprintf(buffer_in, "%X,\0", sigCmpx[n].real);
putsUART1((unsigned int*)buffer_in);//送信する文字列をUART送信バッファに書き込む関数。NULL文字で終了する。
sprintf(buffer_in, "%X,\0\n", sigCmpx[n].imag);
putsUART1((unsigned int*)buffer_in);//送信する文字列をUART送信バッファに書き込む関数。NULL文字で終了する。
}
strcpy(buffer_in,"--fft--,\0");
putsUART1((unsigned int*)buffer_in);//送信する文字列をUART送信バッファに書き込む関数。NULL文字で終了する。
//twiddleFactorsは定義済みtwiddleFactors.c
FFTComplexIP (LOG2_BLOCK_LENGTH, &sigCmpx[0], (fractcomplex *) __builtin_psvoffset(&twiddleFactors[0]), (int) __builtin_psvpage(&twiddleFactors[0]));
/* Store output samples in bit-reversed order of their addresses */
BitReverseComplex (LOG2_BLOCK_LENGTH, &sigCmpx[0]);
for (n=0 ; n < FFT_BLOCK_LENGTH ;n++)
{
sprintf(buffer_in, "%X,\0", sigCmpx[n].real);
putsUART1((unsigned int*)buffer_in);//送信する文字列をUART送信バッファに書き込む関数。NULL文字で終了する。
sprintf(buffer_in, "%X,\0", sigCmpx[n].imag);
putsUART1((unsigned int*)buffer_in);//送信する文字列をUART送信バッファに書き込む関数。NULL文字で終了する。
}
while (1)
{}
return 0;
}
// ******************************************************************************************************
//*******************************************************************************************************
参考文献
(1)トランジスタ技術2007年8月号,9月号
(2)「DSP関数を使ってみよう」トランジスタ技術2008年3月号,4月号,5月号
(3)「Excelを使ったOFDMのシミュレーション」トランジスタ技術2008年11月号pp153-157
(4)DSPライブラリでFIRフィルターを
(5)FFT bug/documentation error