Hi everyone, I am currently doing my Final Year Project. I have coded the
algorithm used for Simulation. It is based on simple Subband Adaptive Filtering
with 2 Bands. One of the Analysis Banks contains low pass filter as its
coefficient and the other contains high Pass filter as its coefficient. (I will
post up the diagram if it is helpful).
This part is under process section
Immediately after the arrays (x0, x1, out_sys0, out_sys1) come out of the
decimation function, everything seems fine. However, these data are gone after
passing into the 'for' loop. Is there anyway i am able to retain the
data? I tried to store the data into another array immediately after the
decimation function but invalid.
I tried to remove the 'for' loop and place the interpolation function
right after the decimation function. The program is working.
In addition, I used fir1 function in Matlab to generate the coefficient for
analysis banks and synthesis banks. However, I am unable to generate even
numbers for my high pass filter i.e. fir1(127,0.5,'high'). Is there
any way I can generate even coefficient for high pass filter? I need this
because for Interpolation,
fir_init(interp_state1, coef_interp1, delay_interp1, 128/2, 2)
I need the coefficient to be even else it will give me decimal.
I hope someone out there can help me. My email/MSN is k...@gmail.com. If anyone
has any sample codes on decimation or interpolation, I would gladly have it.
Thanks a lot.
Below contains the content for ale_demo and process:
========================================================================ale_demo
--------
#include
#include
#include
#include
#define NUMPOINTS 1024
#define TAPS 32
#define FIRTAPS 32
section ("Buffer") fract16 input_arr[NUMPOINTS];
section ("Buffer") fract16 out1[NUMPOINTS];
section ("Buffer") fract16 out2[NUMPOINTS];
//Unknown system
fract16 coefFIR[FIRTAPS] {49,56, 13,-95, -200, -134, 190, 570, 554, -164, -1252,
-1694, -413, 2772, 6708, 9424, 9424, 6708, 2772, -413, -1694, -1252, -164, 554,
570, 190, -134, -200, -95, 13, 56, 49};
fir_state_fr16 sys_state;
fract16 firdelay[FIRTAPS]={0};
//Analysis Bank
// A0
fir_state_fr16 dec_state0;
fract16 delay_dec0[129]; //128: follow Text
fract16 coef_dec0[129] {0, -13, 0, 15, 0, -17, 0, 20, 0, -24, 0, 29, 0, -35, 0,
42, 0, -51, 0, 62, 0, -74, 0, 87, 0, -103, 0, 121, 0, -141, 0, 164, 0, -189, 0,
219, 0, -252, 0, 290, 0, -334, 0, 386, 0, -447, 0, 521, 0, -613, 0, 730, 0,
-887, 0, 1108, 0, -1451, 0, 2058, 0, -3461, 0, 10429, 16390, 10429, 0, -3461, 0,
2058, 0, -1451, 0, 1108, 0, -887, 0, 730, 0, -613, 0, 521, 0, -447, 0, 386, 0,
-334, 0, 290, 0, -252, 0, 219, 0, -189, 0, 164, 0, -141, 0, 121, 0, -103, 0, 87,
0, -74, 0, 62, 0, -51, 0, 42, 0, -35, 0, 29, 0, -24, 0, 20, 0, -17, 0, 15, 0,
-13, 0};
// A1
fir_state_fr16 dec_state1;
fract16 delay_dec1[129]; //128: follow Text
fract16 coef_dec1[129] {0, -13, 0, 15, 0, -17, 0, 20, 0, -24, 0, 29, 0, -35, 0,
42, 0, -51, 0, 62, 0, -74, 0, 87, 0, -103, 0, 121, 0, -141, 0, 164, 0, -189, 0,
219, 0, -252, 0, 290, 0, -334, 0, 386, 0, -447, 0, 521, 0, -613, 0, 730, 0,
-887, 0, 1108, 0, -1451, 0, 2058, 0, -3461, 0, 10429, 16390, 10429, 0, -3461, 0,
2058, 0, -1451, 0, 1108, 0, -887, 0, 730, 0, -613, 0, 521, 0, -447, 0, 386, 0,
-334, 0, 290, 0, -252, 0, 219, 0, -189, 0, 164, 0, -141, 0, 121, 0, -103, 0, 87,
0, -74, 0, 62, 0, -51, 0, 42, 0, -35, 0, 29, 0, -24, 0, 20, 0, -17, 0, 15, 0,
-13, 0};
//Systhesis Bank
// B0
fir_state_fr16 interp_state0;
fract16 delay_interp0[129]; //128: follow Text
fract16 coef_interp0[129] {0, 13, 0, -15, 0, 17, 0, -20, 0, 24, 0, -29, 0, 35,
0, -42, 0, 51, 0, -62, 0, 74, 0, -87, 0, 103, 0, -121, 0, 141, 0, -164, 0, 189,
0, -219, 0, 252, 0, -290, 0, 334, 0, -386, 0, 447, 0, -521, 0, 613, 0, -730, 0,
887, 0, -1108, 0, 1451, 0, -2058, 0, 3461, 0, -10429, 16390, -10429, 0, 3461, 0,
-2058, 0, 1451, 0, -1108, 0, 887, 0, -730, 0, 613, 0, -521, 0, 447, 0, -386, 0,
334, 0, -290, 0, 252, 0, -219, 0, 189, 0, -164, 0, 141, 0, -121, 0, 103, 0, -87,
0, 74, 0, -62, 0, 51, 0, -42, 0, 35, 0, -29, 0, 24, 0, -20, 0, 17, 0, -15, 0,
13, 0};
// B1
fir_state_fr16 interp_state1;
fract16 delay_interp1[129]; //128: follow Text
fract16 coef_interp1[129] {0, -13, 0, 15, 0, -17, 0, 20, 0, -24, 0, 29, 0, -35,
0, 42, 0, -51, 0, 62, 0, -74, 0, 87, 0, -103, 0, 121, 0, -141, 0, 164, 0, -189,
0, 219, 0, -252, 0, 290, 0, -334, 0, 386, 0, -447, 0, 521, 0, -613, 0, 730, 0,
-887, 0, 1108, 0, -1451, 0, 2058, 0, -3461, 0, 10429, 16390, 10429, 0, -3461, 0,
2058, 0, -1451, 0, 1108, 0, -887, 0, 730, 0, -613, 0, 521, 0, -447, 0, 386, 0,
-334, 0, 290, 0, -252, 0, 219, 0, -189, 0, 164, 0, -141, 0, 121, 0, -103, 0, 87,
0, -74, 0, 62, 0, -51, 0, 42, 0, -35, 0, 29, 0, -24, 0, 20, 0, -17, 0, 15, 0,
-13, 0};
//Filtering
//W0
fract16 coefLMS0[TAPS/2];
fract16 lmsdelay0[TAPS/2+NUMPOINTS-1];
fir_state_fr16 firlmsstate0;
fract16 firlmsdelay0[TAPS/2];
//W1
fract16 coefLMS1[TAPS/2];
fract16 lmsdelay1[TAPS/2+NUMPOINTS-1];
fir_state_fr16 firlmsstate1;
fract16 firlmsdelay1[TAPS/2];
////////////////////////
// function prototypes
////////////////////////
void create_sine(fract16*, float);
void add_random(fract16*);
// LMS filter function prototype
void LMS_filter(fract16[], fract16[], fract16[], fract16[], fract16[],
fract16[], fract16[],
int, int, fir_state_fr16*, fir_state_fr16*, fir_state_fr16*,
fir_state_fr16*, fir_state_fr16*,
fir_state_fr16*, fir_state_fr16*);
void main()
{
float freq = 3000;
int i;
//Initializing the unknown system
fir_init(sys_state, coefFIR, firdelay, FIRTAPS, 0);
//Initializing the Analysis Bank
fir_init(dec_state0, coef_dec0, delay_dec0, 128, 2); //fir_init_dec0 - 128:
Follow TEXT
fir_init(dec_state1, coef_dec1, delay_dec1, 128, 2); //fir_init_dec1
//Initializing the Systhesis Bank
fir_init(interp_state0, coef_interp0, delay_interp0, 128/2, 2);
//fir_init_dec0 - 128: Follow TEXT
fir_init(interp_state1, coef_interp1, delay_interp1, 128/2, 2);
//fir_init_dec1
//W0 , W1
fir_init(firlmsstate0, coefLMS0, firlmsdelay0, TAPS/2, 0);
fir_init(firlmsstate1, coefLMS1, firlmsdelay1, TAPS/2, 0);
for (i=0; i
{
create_sine(&input_arr[i], freq);
// add_random(&input_arr[i]);
}
// do processing
LMS_filter(input_arr,out1, out2, lmsdelay0, lmsdelay1, coefLMS0, coefLMS1,
TAPS,NUMPOINTS,&sys_state, &dec_state0, &dec_state1, &interp_state0,
&interp_state1, &firlmsstate0, &firlmsstate1);
printf("Completion of ALE \n");
//}
}
void create_sine(fract16* out, float f)
{
static int j = 0;
float fs = 48000.0;
float step;
step = (float)j/fs;
*out = (fract16) (0.5*sin(2*PI*f*step)*32768.0);
j = __builtin_circindex(j, 1, 48000);
}
void add_random(fract16* out)
{
// As rand() gives value in the range [0, 2^30-1],
// proper scaling and offsetting are implemented.
*out += (fract16)((rand() - 0x20000000)>>16);
}
========================================================================Process
-------
#ifndef ETSI_SOURCE
#define ETSI_SOURCE
#endif
#include
#include
#include
#include
#define BETA 0x0400 //0x0080 // 1/256
#define BETA1 0x7c00 //0x7F80 // 1-(1/256)
#define NUM 1024
fract16 Calc_mu(fract16[], int);
void LMS_filter(fract16 in[],fract16 out1[],fract16 out2[],fract16 lmsdelay0[],
fract16 lmsdelay1[], fract16 coefLMS0[], fract16 coefLMS1[],
int taps,int size, fir_state_fr16* sys_state,
fir_state_fr16* dec_state0, fir_state_fr16* dec_state1,
fir_state_fr16* interp_state0, fir_state_fr16* interp_state1,
fir_state_fr16* firlmsstate0, fir_state_fr16* firlmsstate1)
{
fract16 acc = 0, err = 0, fir_in = 0, mu, mu_err;
int j;
fract16 sys = 0, y0, y1, err0, err1, mu0, mu1, mu_err0, mu_err1;
int i;
int subs = size/2-1;
fract16 sys_array[NUM] = {0};
fract16 out_sys0[NUM/2] = {0};
fract16 out_sys1[NUM/2] = {0};
fract16 out_sys_bank1[NUM/2] = {0};
fract16 x0[NUM/2] = {0};
fract16 x1[NUM/2] = {0};
fract16 err0_array[NUM/2] = {0};
fract16 err1_array[NUM/2] = {0};
fract16 err0_array_out[NUM] = {0};
fract16 err1_array_out[NUM] = {0};
for (i=0; i
{
fir_fr16(&in[i], &sys, 1, sys_state); //fir_fr16(&in[i], &fir_in, 1, s);
sys_array[i] = sys;
// out1[i] = sys_array[i];
}
//Unknown system
fir_decima_fr16(sys_array, out_sys0, size, dec_state0); //dec0
fir_decima_fr16(sys_array, out_sys1, size, dec_state1); //dec1
//Adaptive Filtering - Analysis Bank
fir_decima_fr16(in, x0, size, dec_state0);
fir_decima_fr16(in, x1, size, dec_state1);
for (i=0; i
{
lmsdelay0[size/2-i-1] = x0[i];
lmsdelay1[size/2-i-1] = x1[i];
// call fir filter built-in function
fir_fr16(&lmsdelay0[subs-i], &y0, 1, firlmsstate0);
fir_fr16(&lmsdelay1[subs-i], &y1, 1, firlmsstate1);
err0 = out_sys0[i]-y0;
err1 = out_sys1[i]-y1;
err0_array[i] = err0;
err1_array[i] = err1;
mu0 = Calc_mu(&lmsdelay0[subs-i], taps); // compute step size based on power
of input signal
mu1 = Calc_mu(&lmsdelay1[subs-i], taps); // compute step size based on power
of input signal
//mu = 0x800; //0x01, 0x10, 0x20, 0x100
mu_err0 = multr_fr1x16(mu0, err0);
mu_err1 = multr_fr1x16(mu1, err0);
for (j=0; j
{
coefLMS0[j] = coefLMS0[j] + multr_fr1x16(mu_err0,
multr_fr1x16(lmsdelay0[(subs-i)+j],lmsdelay0[(subs-i)+j]));
coefLMS1[j] = coefLMS1[j] + multr_fr1x16(mu_err1,
multr_fr1x16(lmsdelay1[(subs-i)+j],lmsdelay1[(subs-i)+j]));
}
}
//Adaptive Filtering - Systhesis Bank
fir_interp_fr16(err0_array,err0_array_out,size/2, interp_state0);
fir_interp_fr16(err1_array,err1_array_out,size/2, interp_state1);
for (i=0; i
{
out1[i] = lmsdelay0[i];
}
}
fract16 Calc_mu(fract16 current[], int n)
{
static fract16 est_power = 0; // also used to store previous estimation of
power
// P(n) = (1-Beta) * P(n-1) + Beta * (x(n)^2)
est_power = multr_fr1x16
(
BETA1,
est_power
)
+
multr_fr1x16
(
BETA,
multr_fr1x16
(
current[0],
current[0]
)
);
// The following computes: (0.25/NPx)
// Precomputed constant: 0.25/32 (or 0x100 in fract16)
return div_s(0x800, est_power);
}
========================================================================
Some problems with Decimation/Interpolation for Subband Adaptive Filtering (Adaptive Filtering)
Started by ●April 13, 2009
