diff --git a/test/src/IIR_Filter/.vscode/settings.json b/test/src/IIR_Filter/.vscode/settings.json new file mode 100644 index 0000000..4727727 --- /dev/null +++ b/test/src/IIR_Filter/.vscode/settings.json @@ -0,0 +1,6 @@ +{ + "files.associations": { + "filter.h": "c", + "coefs.h": "c" + } +} \ No newline at end of file diff --git a/test/src/IIR_Filter/400Hz_bandpass.fda b/test/src/IIR_Filter/400Hz_bandpass.fda new file mode 100644 index 0000000..0809669 Binary files /dev/null and b/test/src/IIR_Filter/400Hz_bandpass.fda differ diff --git a/test/src/IIR_Filter/CMakeLists.txt b/test/src/IIR_Filter/CMakeLists.txt new file mode 100644 index 0000000..41d08ab --- /dev/null +++ b/test/src/IIR_Filter/CMakeLists.txt @@ -0,0 +1,5 @@ +cmake_minimum_required(VERSION 3.12) +project(iirfilter) +aux_source_directory(. SOURCE) +message("source " ${SOURCE}) +add_executable(iirfilter ${SOURCE}) diff --git a/test/src/IIR_Filter/Makefile b/test/src/IIR_Filter/Makefile new file mode 100644 index 0000000..f810be0 --- /dev/null +++ b/test/src/IIR_Filter/Makefile @@ -0,0 +1,23 @@ +TARGET = filter + +HEADERS = $(wildcard *.h) +SRC = $(wildcard *.c) +OBJS = $(SRC:.c=.o) + +CC = gcc +#CFLAGS = -Wall -std=gnu99 + +#can be deleted, .c can be translated to .o automatically +#use this to "@echo" +%.o:%.c + @$(CC) -g -o $@ -c $< + +$(TARGET) : $(OBJS) + @$(CC) -g -o $@ $^ + +$(OBJS) : $(HEADERS) + +.PHONY : clean +clean : + rm $(TARGET) $(OBJS) + diff --git a/test/src/IIR_Filter/README.md b/test/src/IIR_Filter/README.md new file mode 100644 index 0000000..9c2c679 --- /dev/null +++ b/test/src/IIR_Filter/README.md @@ -0,0 +1,64 @@ +## An IIR filter designed with MATLAB and implemented in C + +More details can be found in [my blog](https://www.cnblogs.com/lyrich/p/10987875.html). + +The transfer function of a biquad filter can be written as: H(z)=(b₀+b₁z⁻¹+b₂z⁻²)/(a₀+a₁z⁻¹+a₂z⁻²). + +Use Filter Designer in MATLAB to design a 400 Hz bandpass IIR filter: +![在这里插入图片描述](https://img-blog.csdnimg.cn/20190601212647210.png) + +Analysis->Filter Coeffients : +![在这里插入图片描述](https://img-blog.csdnimg.cn/20190601213027484.png) +Numerator,from top to bottom: b₀,b₁ and b₂. +Denominator, from top to bottom: a₀,a₁ and a₂,a₀ is always 1. + +Use arrays to store filter coeffients: +```c +//8-order IIR filter with 4 sections +const int sections = 4; + +//nominator +const float b[4][3] = { + { 1, -1.984454421, 1 }, + { 1, -1.999405318, 1 }, + { 1, -1.993167556, 1 }, + { 1, -1.998644244, 1 } +}; + +//denominator +const float a[4][3] = { + { 1, -1.984532740, 0.9884094716 }, + { 1, -1.988571923, 0.9909378613 }, + { 1, -1.991214225, 0.9962624248 }, + { 1, -1.995917854, 0.9977478940 } +}; + +const float gain[4] = { 0.583805364, 0.583805364, 0.170388576, 0.170388576 }; +``` + +Store filter states: +```c +float w[sections][2]; //filter states +``` + +Initialize filter: +```c +for (int i = 0; i < sections; i++) { + w[i][0] = 0; //w[n-1] + w[i][1] = 0; //w[n-2] +} +``` + +Calculate filter output: +```c +y[0] = pcmIn[i]; + for (j = 0; j < sections; j++) { + tmp[j] = y[j] - a[j][1] * w[j][0] - a[j][2] * w[j][1]; //calculate w[n] + y[j+1] = tmp[j] + b[j][1] * w[j][0] + b[j][2] * w[j][1]; //calculate the j-th section filter output y[n] + w[j][1] = w[j][0]; //move w[n-1] -> w[n-2] + w[j][0] = tmp[j]; //move w[n] -> w[n-1] + y[j+1] = gain[j] * y[j+1]; //multiply with gain + } + +out = y[j]; +``` diff --git a/test/src/IIR_Filter/coefs.c b/test/src/IIR_Filter/coefs.c new file mode 100644 index 0000000..d37beb9 --- /dev/null +++ b/test/src/IIR_Filter/coefs.c @@ -0,0 +1,107 @@ +#include "coefs.h" + +const int sections = 4; + +const float b[4][3] = { + { 1, -1.984454421, 1 }, + { 1, -1.999405318, 1 }, + { 1, -1.993167556, 1 }, + { 1, -1.998644244, 1 } +}; + +//denominator +const float a[4][3] = { + { 1, -1.984532740, 0.9884094716 }, + { 1, -1.988571923, 0.9909378613 }, + { 1, -1.991214225, 0.9962624248 }, + { 1, -1.995917854, 0.9977478940 } +}; + +const float gain[4] = { 0.583805364, 0.583805364, 0.170388576, 0.170388576 }; + + +#define MWSPT_NSEC 13 +const int NL[MWSPT_NSEC] = { 1,3,1,3,1,3,1,3,1,3,1,3,1 }; +const float NUM[MWSPT_NSEC][3] = { + { + 0.01416839943433, 0, 0 + }, + { + 1, 0, -1 + }, + { + 0.01416839943433, 0, 0 + }, + { + 1, 0, -1 + }, + { + 0.01061464435972, 0, 0 + }, + { + 1, 0, -1 + }, + { + 0.01061464435972, 0, 0 + }, + { + 1, 0, -1 + }, + { + 0.005015094751417, 0, 0 + }, + { + 1, 0, -1 + }, + { + 0.005015094751417, 0, 0 + }, + { + 1, 0, -1 + }, + { + 0.8912509381337, 0, 0 + } +}; +const int DL[MWSPT_NSEC] = { 1,3,1,3,1,3,1,3,1,3,1,3,1 }; +const float DEN[MWSPT_NSEC][3] = { + { + 1, 0, 0 + }, + { + 1, -1.992740305702, 0.9977910469633 + }, + { + 1, 0, 0 + }, + { + 1, -1.996837669615, 0.9986690342012 + }, + { + 1, 0, 0 + }, + { + 1, -1.989872067756, 0.9942825585574 + }, + { + 1, 0, 0 + }, + { + 1, -1.993969383703, 0.9960603774853 + }, + { + 1, 0, 0 + }, + { + 1, -1.989477089804, 0.9929580678414 + }, + { + 1, 0, 0 + }, + { + 1, -1.991214772896, 0.9938596848053 + }, + { + 1, 0, 0 + } +}; diff --git a/test/src/IIR_Filter/coefs.h b/test/src/IIR_Filter/coefs.h new file mode 100644 index 0000000..2d79553 --- /dev/null +++ b/test/src/IIR_Filter/coefs.h @@ -0,0 +1,15 @@ +#ifndef _COEFS_H +#define _COEFS_H + +//8-order IIR filter with 4 sections +extern const int sections; + +//nominator +extern const float b[4][3]; + +//denominator +const float a[4][3]; + +const float gain[4]; + +#endif diff --git a/test/src/IIR_Filter/filter.c b/test/src/IIR_Filter/filter.c new file mode 100644 index 0000000..116f696 --- /dev/null +++ b/test/src/IIR_Filter/filter.c @@ -0,0 +1,38 @@ +#include "coefs.h" +#include "filter.h" + +void filter_reset(float w[][2]) +{ + for (int i = 0; i < sections; i++) { + w[i][0] = 0; + w[i][1] = 0; + } +} + +void filter(short* pcmIn, short* pcmOut, float w[][2], int sample_size) +{ + float y[sections+1], tmp[sections], out; + int i,j; + + for (i = 0; i < sample_size; i++) { + y[0] = pcmIn[i]; + for (j = 0; j < sections; j++) { + tmp[j] = y[j] - a[j][1] * w[j][0] - a[j][2] * w[j][1]; //calculate w[n] + y[j+1] = tmp[j] + b[j][1] * w[j][0] + b[j][2] * w[j][1]; //calculate the j-th section filter output y[n] + w[j][1] = w[j][0]; //move w[n-1] -> w[n-2] + w[j][0] = tmp[j]; //move w[n] -> w[n-1] + y[j+1] = gain[j] * y[j+1]; //multiply with gain + } + + out = y[j]; + + if (out > 32767) + out = 32767; + + if (out < -32768) + out = -32768; + + pcmOut[i] = (short)(out); + } +} + diff --git a/test/src/IIR_Filter/filter.h b/test/src/IIR_Filter/filter.h new file mode 100644 index 0000000..a55e4a2 --- /dev/null +++ b/test/src/IIR_Filter/filter.h @@ -0,0 +1,25 @@ +#ifndef _FILTER_H +#define _FILTER_H + +//WAVE PCM soundfile format +typedef struct +{ + char chunk_id[4]; + int chunk_size; + char format[4]; + char subchunk1_id[4]; + int subchunk1_size; + short int audio_format; + short int num_channels; + int sample_rate; //sample_rate denotes the sampling rate. + int byte_rate; + short int block_align; + short int bits_per_sample; + char subchunk2_id[4]; + int subchunk2_size; //subchunk2_size denotes the number of samples. +} wav_header; + +void filter_reset(float w[][2]); +void filter(short* pcmIn, short* pcmOut, float w[][2], int sample_size); + +#endif diff --git a/test/src/IIR_Filter/in.wav b/test/src/IIR_Filter/in.wav new file mode 100644 index 0000000..f5c5812 Binary files /dev/null and b/test/src/IIR_Filter/in.wav differ diff --git a/test/src/IIR_Filter/main.c b/test/src/IIR_Filter/main.c new file mode 100644 index 0000000..e3bdee1 --- /dev/null +++ b/test/src/IIR_Filter/main.c @@ -0,0 +1,36 @@ +#include +#include "filter.h" +#include "coefs.h" +#define WSIZE 512 + +short bufferIn[WSIZE] = {0}; +short bufferOut[WSIZE] = {0}; + +int main() +{ + FILE * fp_in, *fp_out; + int i = 0, nb; + wav_header wav_header_buf; + + float w[sections][2]; //filter states + + fp_in = fopen("in.wav", "rb"); + fp_out = fopen("out.wav", "wb"); + + fread(&wav_header_buf, 1, sizeof(wav_header), fp_in); + fwrite(&wav_header_buf,1, sizeof(wav_header), fp_out); + + while (!feof(fp_in)) { + + nb = fread(bufferIn, 2, WSIZE, fp_in); + + filter_reset(w); + filter(bufferIn, bufferOut, w, WSIZE); + + fwrite(bufferOut, 2, nb, fp_out); //writing read data into output file + } + + printf("completed\n"); + + return 0; +} diff --git a/test/src/IIR_Filter/out.wav b/test/src/IIR_Filter/out.wav new file mode 100644 index 0000000..1df72a0 Binary files /dev/null and b/test/src/IIR_Filter/out.wav differ