Skip to content

Commit

Permalink
Chebyshev coefficients validated. X and Y sequences stored correctly.…
Browse files Browse the repository at this point in the history 
… Validated by spectrum analyze in python.
  • Loading branch information
Nicolas Gravillon committed Jul 26, 2023
1 parent 850500b commit b0ad467
Show file tree
Hide file tree
Showing 6 changed files with 62 additions and 1,051 deletions.
2 changes: 1 addition & 1 deletion Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,4 +15,4 @@ filter.o:
$(CC) $(CFLAGS) -c filter.c

clean:
\rm *.o example
\rm *.o lpbutter_wav
Binary file added data/guitar_E_filtered.wav
View file
Binary file not shown.
28 changes: 16 additions & 12 deletions filter.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ ChebFilter* call_205(int P, ChebFilter* filter, double FC, int NP, int LH, doubl
double rp= -cos(M_PI/(NP*2) + (P-1)*M_PI/NP);
double ip= sin(M_PI/(NP*2) + (P-1)*M_PI/NP);

if(DEBUG) {
if(DEBUG>1) {
printf("\n[call_205 #%d] rp= %.10lf\n", P, rp);
printf("[call_205 #%d] ip= %.10lf\n", P, ip);
printf("\n");
Expand All @@ -35,7 +35,7 @@ ChebFilter* call_205(int P, ChebFilter* filter, double FC, int NP, int LH, doubl
kx= (exp(kx) + exp(-kx))/2;
rp= rp * ((exp(vx) - exp(-vx))/2.0)/kx;
ip= ip * ((exp(vx) + exp(-vx))/2.0)/kx;
if(DEBUG) {
if(DEBUG>1) {
printf("[call_205 #%d PR!=0] rp= %.10lf\n", P, rp);
printf("[call_205 #%d PR!=0] ip= %.10lf\n", P, ip);
printf("[call_205 #%d PR!=0] es= %.10lf\n", P, es);
Expand All @@ -54,7 +54,7 @@ ChebFilter* call_205(int P, ChebFilter* filter, double FC, int NP, int LH, doubl
double y1=(8 - 2*m*pow(t,2))/d;
double y2=(-4 - 4*rp*t - m*pow(t,2))/d;

if(DEBUG) {
if(DEBUG>1) {
printf("[call_205 #%d] t= %.10lf\n", P, t);
printf("[call_205 #%d] w= %.10lf\n", P, w);
printf("[call_205 #%d] m= %.10lf\n", P, m);
Expand All @@ -74,7 +74,7 @@ ChebFilter* call_205(int P, ChebFilter* filter, double FC, int NP, int LH, doubl
if(LH==0) {
k= sin((double)1/2 - w/2)/sin((double)1/2 + w/2);
}
if(DEBUG) printf("[call_205 #%d] k= %.10lf\n", P, k);
if(DEBUG>1) printf("[call_205 #%d] k= %.10lf\n", P, k);

d= 1 + y1*k - y2*pow(k, 2);
filter->a0= (x0 - x1*k + x2*pow(k,2))/d;
Expand All @@ -83,7 +83,7 @@ ChebFilter* call_205(int P, ChebFilter* filter, double FC, int NP, int LH, doubl
filter->b1= (2*k + y1 + y1*pow(k,2) - 2*y2*k)/d;
filter->b2= (-pow(k,2) - y1*k +y2)/d;

if(DEBUG) {
if(DEBUG>1) {
printf("[call_205 #%d] d= %.10lf\n", P, d);
printf("[call_205 #%d] a0= %.10lf\n", P, filter->a0);
printf("[call_205 #%d] a1= %.10lf\n", P, filter->a1);
Expand All @@ -96,7 +96,7 @@ ChebFilter* call_205(int P, ChebFilter* filter, double FC, int NP, int LH, doubl
filter->b1=-filter->b1;
}

if(DEBUG) {
if(DEBUG>1) {
printf("[call_205 #%d] a1= %.10lf\n", P, filter->a1);
printf("[call_205 #%d] b1= %.10lf\n", P, filter->b1);
}
Expand Down Expand Up @@ -181,7 +181,7 @@ ChebFilter* create_che_filter(int NP, double PR, int LH, double FC) {
filter->b[i]=b[i];
}

if(DEBUG) {
if(DEBUG>0) {
printf("Order: %d\n", NP);
printf("Fc: %lf\n", FC);
printf("a0: %lf\n", filter->a0);
Expand All @@ -191,14 +191,16 @@ ChebFilter* create_che_filter(int NP, double PR, int LH, double FC) {
printf("b2: %lf\n", filter->b2);

printf("a= [");
for(int i=0; i<20; i++) {
printf("%.10e ", a[i]);
for(int i=0; i<=NP; i++) {
//printf("%.10e ", a[i]);
printf("%.10lf ", a[i]);
}
printf("]\n");

printf("b= [");
for(int i=0; i<20; i++) {
printf("%.10e ", b[i]);
for(int i=0; i<=NP; i++) {
//printf("%.10e ", b[i]);
printf("%.10lf ", b[i]);
}
printf("]\n");
}
Expand Down Expand Up @@ -232,7 +234,7 @@ double applyfilter(ChebFilter* filter, double X0) {
filter->X[0]=X0;
double Y=0;
for(int i=0; i<=filter->NP; i++)
Y= filter->a[i]*filter->X[i] + filter->b[i]*filter->Y[i];
Y+= filter->a[i]*filter->X[i] + filter->b[i]*filter->Y[i];

filter->Y[0]=Y;
if(DEBUG) {
Expand All @@ -241,6 +243,8 @@ double applyfilter(ChebFilter* filter, double X0) {
printf("X[%d]: %.10lf\n", i, filter->X[i]);
for(int i=0; i<=filter->NP; i++)
printf("Y[%d]: %.10lf\n", i, filter->Y[i]);

printf("X: %.6lf - Y: %.6lf\n", X0, Y);
}

for(int i=filter->NP; i>0; i--) {
Expand Down
Binary file removed lpbutter_wav
View file
Binary file not shown.
83 changes: 45 additions & 38 deletions lpbutter_wav.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,9 @@

#include "filter.h"

#define DEBUG 1
#define DEBUG 0
#define DATA_LEN 1006507 // # frames
#define RUNS 3
#define RUNS 1
#define ORDER 2


Expand All @@ -16,25 +16,30 @@
#define SAMPLE_RATE 48000
#define BLOCK_SIZE 480

#define FC 205
//#define FC 205
#define FC 102.5


int main() {

//ChebFilter* filter = create_bw_low_pass_filter(ORDER, (double)2*M_PI*205/44100);
/*
NP : Number of poles / Order
PR: Percent ripple
LH: 0 for lowpass, 1 for highpass
FC: cutoff frequency (0 to 0.5), typically fc(Hz)/SR(Hz), with fc<SR/2. Example: (double)1000/44100
*/
//ChebFilter* filter = create_bw_low_pass_filter(ORDER, (double)FC/SAMPLE_RATE);
//ChebFilter* filter = create_bw_low_pass_filter(ORDER, 0.01);
ChebFilter* filter = create_che_filter(ORDER, 0.5, 0, 0.01);
ChebFilter* filter = create_che_filter(4, 0.5, 0, FC/48000.0);
double input[DATA_LEN];
//double filtered_signal[RUNS][DATA_LEN];
double *output;

output=(double*)malloc(DATA_LEN*sizeof(double));

TinyWav tw;
tinywav_open_read(&tw, "data/guitar_E.wav", TW_SPLIT );

//printf("# occurences: %d\n", (int)floor(DATA_LEN/BLOCK_SIZE));

int line=0;
//for (int i = 0; i < (int)floor((double)DATA_LEN/BLOCK_SIZE); i++) {
for (int i = 0; i < 2096; i++) {
float samples[NUM_CHANNELS * BLOCK_SIZE];

Expand All @@ -46,49 +51,51 @@ int main() {

tinywav_read_f(&tw, samplePtrs, BLOCK_SIZE);

for(int k=0; k<BLOCK_SIZE; k++) {
//printf("#%d ", line);
input[line++]=(double)samplePtrs[0][i];
//printf("%0.6f\n", samplePtrs[0][k]);
}
for(int k=0; k<BLOCK_SIZE; k++)
input[line++]=(double)samplePtrs[0][k];
}

tinywav_close_read(&tw);

/*
printf("Opening and reading data file...\n");
FILE *fp = fopen("data/inputs.raw", "r");
if(fp == NULL) {
perror("Unable to open file!");
return(-1);
}
char line[30];
unsigned int i=0;
while(fgets(line, sizeof(line), fp) != NULL) {
sscanf(line, "%lf", &input[i]);
i++;
}
fclose(fp);
printf("Data stored into input[].\n");
*/

for(int i=0; i<line; i++)
output[i] = applyfilter(filter, input[i]);

/*
for(int j=0; j<line; j++) {
printf("Run #%d:\n", j);
for(int i=0; i<DATA_LEN; i++) {
filtered_signal[j][i] = applyfilter(filter, input[i]);
}
tinywav_open_write(&tw, NUM_CHANNELS, SAMPLE_RATE,
TW_FLOAT32, // the output samples will be 32-bit floats. TW_INT16 is also supported
TW_INLINE, // the samples to be written will be assumed to be inlined in a single buffer.
// Other options include TW_INTERLEAVED and TW_SPLIT
"data/guitar_E_filtered.wav" // the output path
);
int c=0;
for (int i = 0; i < 2096; i++) {
float samples[480 * NUM_CHANNELS];
for(int i=0; i<480; i++)
samples[i]=(float)output[c++];
tinywav_write_f(&tw, samples, sizeof(samples));
}
tinywav_close_write(&tw);
*/
for(int i=0; i<line; i++)
printf("%.10lf\n", output[i]);


/*
printf("Comparing results:\n");
for(int i=0; i<500; i++) {
printf("[bp_filter %i] input: %lf | ", i, input[i]);
for(int j=0; j<RUNS; j++)
printf(" f_out#%d %lf |", j, filtered_signal[j][i]);
printf(" f_out#1 %lf |", output[i]);
printf(" a= [");
for(int j=0; j<3; j++)
printf(" %.6lf ", filter->a[j]);
printf("] | b= [");
for(int j=0; j<3; j++)
printf(" %.6lf ", filter->b[j]);
printf("]");
printf("\n");
}
*/
Expand Down
Loading

0 comments on commit b0ad467

Please sign in to comment.