📄 real_dit_fft.c
字号:
/*******************************************************************************
; Module: Real_DIT_FFT
; Filename: Real_DIT_FFT.c
; Project: DSP library for XC16x Microcontroller
;------------------------------------------------------------------------------
; Compiler: Keil
;
; Version: V1.2
;
; Description: Implementation of forward real Radix-2 decimation-in-time
; Fast Fourie Transformation
;
; Date: May 2003
;
; Copyright: Infineon Technologies AG
;
; Author: Guangyu Wang
;******************************************************************************/
/******************************************************************************
; void real_DIT_FFT(
; DataS* x, // input vector
; DataS* index, // bit reversed input indecies
; DataS exp, // exponent of vector size
; DataS* table, // trigonomic function table
; DataS* X // output of FFT
; )
;
; INPUTS:
; x the input value in 1Q15 format
; index bit reversed input indecie
; exp exponent of vector size
; table trigonomic function table
;
; RETURN:
; X the FFT output vector in 1Q15 format
;
; ALGORITHM: Radix-2 decimation-in-time FFT
;
; REGISTER USAGE:
; 1. From .c file to .asm file:
; R8 pointer of input vector x
; R9 pointer of index vector
; (R10) = exp
; R11 pointer of trigonomic function table
; R12 contains the pointer of output vector X
;
;
; 2. From .asm file to .c file
;
; Assumptions:
; 1. scale factor = 1/N (N=2^exp), preventing overflow
; 2. FFT_of_x(k) = N*X(k)
;
;*****************************************************************************/
/*========================= Define registers =================================
Outloop LIT 'R4' ; R4 counter
Midloop LIT 'R3' ; Midloop counter
Inloop LIT 'R1' ; Inloop counter
P_off LIT 'QR0' ; offset to P element in x
Q_off LIT 'QR1' ; offset to Q element in x
FFT_in LIT 'R8' ; Pointer to first element in x
index LIT 'R9' ; Pointer to first element in index
exp LIT 'R10' ; exp
table LIT 'R11' ; Pointer to first element in table
X LIT 'R12' ; Pointer to first element in X
counter LIT 'R13' ; 2^(exp)=N
============================================================================*/
/* ---------------------------------------------------------
; Butterfly macros used for radix-2 forward
; decimation-in-time (DIT) FFT
; ---------------------------------------------------------
;**************************************************************
;
; n = stage
;
; Pn = Re(Pn) + Im(Pn)
;
; Qn = Re(Qn) + Im(Qn)
;
; Pn------------------->-------->(+)----------------->Pn+1
; \ ^
; \ /
; \ /
; \/
; /\
; / \
; / \
; / V
; Qn------------------->-------->(+)----------------->Qn+1
; k
; W -1
; N
;
; k
; Pn+1 = Pn + W * Qn
; N
; k
; Qn+1 = Pn - W * Qn
; N
;
; k -j(2pi/N)k
; W = e = cos(X) - j * sin(X)
; N
;
; 2*pi
; X = ( ---- )k
; N
;
; Pn+1 = Re(Pn) + Re(Qn) * cos(x) + Im(Qn) * sin(x)
; + j[Im(Pn) + Im(Qn) * cos(x) - Re(Qn) * sin(x)]
;
; Qn+1 = Re(Pn) - Re(Qn) * cos(x) - Im(Qn) * sin(x)
; + j[Im(Pn) - Im(Qn) * cos(x) + Re(Qn) * sin(x)]
;----------------------------------------------------------*/
#include "DspLib_Keil.h"
void real_DIT_FFT(DataS* x, DataS* index, DataS exp, DataS* table, DataS* X)
{
__asm
{
//store the system state
PUSH R13
PUSH R14
PUSH R15
//initialization of registers
MOV R4,exp // (R4) = exp
SUB R4,#1h // (R4) = exp-1
MOV R13,#1h // (R13)=counter = 1
SHL R13,exp // (R13)=counter = 2^exp = N
MOV R5,R13 // (R5) = N
SHR R5,#1h // (R5) = N/2
//initialize the basic address of cosinus
ADD R5,table // (R5) = table + N/2= address of cos(x)
//reset R2
MOV R2,R13 // (R2) = N
SHL R2,#1h // (R2) = 2*N
//MAC registers initialization
MOV MCW,#0200h // MP=0, MS=1
OUT_LOOP:
MOV R3,#0 //(R3) = 0
MOV R6,R13 //(R6) = N
SHR R6,#2 //(R6) = N/4 = number of Inloops
MID_LOOP:
MOV R1,R6 // restore Inloop
EXTR #3
MOV QR0,R3 // (QR0) = offset of Pn in FFT_in
MOV QR1,QR0 // (QR1) = (QR0)
ADD QR1,R13 // (QR1) = offset of Qn in FFT_in
PUSH R13
//------Determination of Index for Twiddle-Factor
MOV R7,R3
SHR R7,R4
//------Determination of Twiddle-Factor
ADD R7,index // (R7)=(R7)+index
MOV R7,[R7] // (R7) = Bit reverse index
MOV R14,R5 // move basic address of cos(x) into R14
MOV R15,table // move basic address of sin(x) into R15
ADD R14,R7 // add offset
ADD R15,R7
MOV R14,[R14] // (R14) = cos(x)
MOV R15,[R15] // (R15) = sin(x)
PUSH R5 // push basic address of cos(x)into stack
PUSH R6
PUSH R12
IN_LOOP:
//------use general butterfly
CoNOP [x+QR0]
// %Wk_Butterfly()
/*-----------------------------------------
; %DEFINE (Wk_Butterfly ())(
; general Butterfly
; (R14) = cos(x)
; (R15) = sin(x)
; W = cos(x) - j * sin(x)
------------------------------------------*/
//------Input Pn---------------------------
//(R5) = Re(Pn)
//(R6) = Im(Pn)
MOV R5,[x+] // (R5) = X[QR0] = Re(Pn)
// (R8) = (R8)+2
MOV R6,[x] // (R6) = X[QR0+1] = Im(Pn)
ASHR R5,#1 // Re(Pn)=Re(Pn)/2,
ASHR R6,#1 // Im(Pn)=Im(Pn)/2,
SUB x,#2h // (R8) = (R8)-2
CoNOP [x - QR0] // (R8) = (R8) - (QR0)
CoNOP [x + QR1] // (R8) = (R8) + (QR1)
//------Input Qn---------------------------
// ((R12)) = Re(Qn)
// ((R12+2)) = Im(Qn)
//------Calculation of Real Parts---------
//(R7) = Re(Pn+1) = Re(Pn) + Re(Qn) * cos(X) + Im(Qn) * sin(x)
MOV MAH,R5 // (ACC) = Re(Pn)
CoMAC R14,[x+] // (ACC) = Re(Pn) + (Re(Qn) * cos(x))
// (R8) = (R8)+2
CoMAC R15,[x-] // (ACC) = (ACC) + (Im(Qn) * sin(x))
// (R8) = (R8) - 2
CoSTORE R7,MAS // (R7) = Re(Pn+1)
//(R8) = Re(Qn+1) = Re(Pn) - Re(Qn) * cos(x) - Im(Qn) * sin(x)
MOV MAH,R5 // (ACC) = Re(Pn)
CoMAC- R14,[x+] // (ACC) = Re(Pn) - (Re(Qn) * cos(x))
// (R8) = (R8) + 2
CoMAC- R15,[x] // (ACC) = (ACC) - (Im(Qn) * sin(x))
CoSTORE R12,MAS // (R12) = Re(Qn+1)
//------Calculation of image parts---------
//(R5) = Im(Pn+1) = Im(Pn) + Im(Qn) * cos(x) - Re(Qn) * sin(x)
MOV MAH,R6 // (ACC) = Im(Pn)
CoMAC R14,[x-] // (ACC) = Im(Pn) + (Im(Qn) * cos(x))
// (R8) = (R8)-2
CoMAC- R15,[x+] // (ACC) = (ACC) - (Re(Qn) * sin(x))
// (R8) = (R8)-2
CoSTORE R5,MAS // (R5) = Im(Pn+1)
//(R9) = Im(Qn+1) = Im(Pn) - Im(Qn) * cos(x) + Re(Qn) * sin(x)
MOV MAH,R6 // (ACC) = Im(Pn)
CoMAC- R14,[x-] // (ACC) = Im(Pn) - (Im(Qn) * cos(x))
// (R8) = (R8)-2
CoMAC R15,[x+] // (ACC) = (ACC) + (Re(Qn) * sin(x))
// (R8) = (R8)+2
CoSTORE R13,MAS // (R13) = Im(Qn+1)
//------Output Pn+1 and Qn+1-----------
MOV [x],R13 // X[QR1+2] = R13 = Im(Qn+1)
MOV [-x],R12 // X[QR1] = R12 = Re(Qn+1)
CoNOP [x - QR1] // (R8) = (R8) - (QR1)
CoNOP [x + QR0] // (R8) = (R8) + (QR0)
ADD x,#2 // (R8) = (R8)+2
MOV [x],R5 // X[QR0+2] = R5 = Im(Pn+1)
MOV [-x],R7 // X[QR0] = R7 = Re(Pn+1)
CoNOP [x - QR0] // (R8) = (R8) - (QR0)
//------Incrementing FFT_in Index Counter--------
EXTR #2
ADD QR0,#4
ADD QR1,#4
//------Decrementing IN_LOOP Counter-------------
SUB R1,#1
JMPR cc_NZ,IN_LOOP
POP R12 // pop basic address of cos(x) out stack
POP R6
POP R5
MOV R13,#1 //(R13)=1
SHL R13,exp //(R13)=N
SHL R13,#1 //(R13)=2*N
//------Increment FFT_IN Index offset
ADD R3,R2
CMP R3,R13
POP R13
JMPR cc_ULT,MID_LOOP
//END_MID_LOOP:
MOV R2,R13 // (R2) = counter
SHR R13,#1h // (R13) = (R13)/2
SUB R4,#1h // R4 = R4 - 1
JMPR cc_NZ,OUT_LOOP
// ***** End Stage 1 - (exp-1) *****
// ********** Last Stage **********
// This stage performs the extraction of the real valued FFT
// out of the previously performed complex FFT. The
// output is sorted in linear-order. Because of symmetric property of
// FFT only the first N/2 points and the Nyquist point are calculated.
//------Extraction of the frequency values with indices [0, N/2]
MOV R2,#(-2) // init UN_LOOP counter
MOV R13,#1h // (R13) = 1
SHL R13,exp // (R13) = 2^(exp) = N
UN_LOOP:
PUSH R5 // push cos(x) into stack
PUSH table
ADD R2,#2 // (R2)=(R2)+2
//------determine input index of R(k) and I(k) in FFT_in
MOV R3,R2 // (R3)=(R2)
ADD R3,index // (R3) = (R2) + index
MOV R14,[R3]
ADD R14,x
//------Input-------------------------------
MOV R3,[R14+] // R3 = FFT_in[R3] = R(k)
MOV R14,[R14] // R14 = FFT_in[R3+2] = I(k)
MOV R6,R3 // R6 = R(k)
MOV R7,R14 // R7 = I(k)
CMP R2,#0
JMPR cc_EQ, DC_output //jump to calculation of DC element
//------determine input index of R(N/2-k) and I(N/2-k) in FFT_in
MOV R7,R13 // (R7) = N
SUB R7,R2 // (R7) = N - R2
ADD R7,index // (R7) = (R7) + index
MOV R7,[R7] // (R7) = [N - index]
ADD R7,x // (R7) = (R7) + FFT_in
//------Input-------------------------------
MOV R6,[R7+] // (R6) = R(N/2-k)
MOV R7,[R7] // (R7) = I(N/2-k)
//------save register-----------------------
DC_output:
MOV R15,R6 // R15 = R6= R(N/2-k)
MOV R1,R14 // R1 = R14= I(k)
//------calculation of the spectra H(k) and G(k)
SUB R15,R3 // (R15) = (R15)-(R3) = -(R(k)-R(N/2-k))=Im{G}
ADD R1,R7 // (R1) = (R1)+(R7) = I(k)+I(N/2-k)=Re{G}
ADD R3,R6 // (R13) = (R13)+(R6) = R(k)+R(k+N/2)=Re{H}
SUB R14,R7 // (R14) = (R14)-(R7) = I(k)-I(N/2-k)=Im{H}
ASHR R15,#1
ASHR R1,#1
ASHR R3,#2 //to make sure same point position in calculation of Re{X(k)}
ASHR R14,#2 //and Im{X(k)}
//------assembly of the full spectrum spectra----
ADD table,R2
ADD R5,R2
MOV R5,[R5] // (R5) = cos(x)
MOV table,[table] // (table) = sin(x)
//------calculate Re{X(k)}=Re{H(k)}+cos(X)*Re{G(k)}+sin(X)*Im{G(k)}
MOV MAH,R3 //(ACC)=(R13)=Re{H(k)}
MOV [R12],R1 //((X))=Re{G}
CoMAC R5,[R12] //(ACC)=(ACC)+cos(x)*Re{G(k)}
MOV [R12],R15 //(X)=(R1) =Im{G}
CoMAC table,[R12] //(ACC)=(ACC)+sin(X)*Im{G(k)}
//Write the real part Re{X(k)} into X
COSTORE [R12+],MAS //((R12))=limited(ACC) =Re{X(k)}
//(R12)=(R12)+2
//-------Calculate Im{X(k)}=Im{H(k)}+cos(x)*Im{G(k)}-sin(x)*Re{G(k)}
MOV MAH,R14 //(ACC)=(R6)=Im{H(k)}
MOV [R12],R15 //((R12))=(R15)=Im{G}
CoMAC R5,[R12] //(ACC)=(ACC)+cos(x)*Im{G(k)}
//(IDX0)=(IDX0)-(QX0)
MOV [R12],R1 //((R12))=Re{G}
CoMAC- table,[R12] //(ACC)=(ACC)-sin(X)*Re{G(k)}
//Write the image part Im{X(k)} into X
COSTORE [R12+],MAS //((R12))=limited(ACC)=Im{X(k)}
//(R12)=(R12)+2
CMP R2,#0
JMPR cc_NE, loop_counter //jump to loop counter
//The folowing instructions computer the FFT value in Nyquist point,
//that is, k=N/2
//------ FFT value in Nyquist popint N/2 ----------------
//------ Re{X(N/2)}=Re{H(0)}-cos(pi*n)*Re{G(0)}-sin(pi*n)*Im{G(0)}
MOV MAH,R3 //(ACC)=(R5)=Re{H(0)}
MOV [R12],R1 //((X))=Re{G}
CoMAC- R5,[R12] //(ACC)=(ACC)-cos(n*pi)*Re{G(0)}
MOV [R12],R15 //((X))=(R7) =Im{G}
CoMAC- table,[R12] //(ACC)=(ACC)-sin(n*pi)*Im{G(0)}
//Write the real part Re{X(N/2)} into X
MOV R6,R13 //(R6)=N
SHL R6,#1 //(R6)=2*N
ADD R6,R12
SUB R6,#4 //(R6)=(X)+2*N-4
COSTORE [R6+],MAS //((R6))=limited(ACC) =Re{X(N/2)}
//(R6)=(R6)+2
//------- Im{X(N/2)}=Im{H(0)}-cos(n*pi)*Im{G(0)}+sin(n*pi)*Re{G(0)}
MOV MAH,R14 //(ACC)=(R6)=Im{H(0)}
MOV [R12],R15 //((X))=(R7)=Im{G(0)}
CoMAC- R5,[R12] //(ACC)=(ACC)-cos(n*pi)*Im{G(0)}
//(IDX0)=(IDX0)-(QX0)
MOV [R12],R1 //((X))=Re{G(0)}
CoMAC table,[R12] //(ACC)=(ACC)+sin(n*pi)*Re{G(0)}
//Write the image part Im{X(k)} into X
COSTORE [R6],MAS //((R6))=limited(ACC)=Im{X(N/2)}
//------ loop counter-------------------------
loop_counter:
POP table
POP R5 // pop the basic address of cos(x)from stack
MOV R1,R13 // (R1)=N
SUB R1,#2 // (R1)=N-2
CMP R2,R1
JMPR cc_ULT,UN_LOOP
//restore the system state
POP R15
POP R14
POP R13
RET
}
}
//------------------- END OF FILE ----------------------------------------------
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -