add iir test

test/src/IIR_Filter/.vscode/settings.json

@@ -0,0 +1,6 @@
+{
+    "files.associations": {
+        "filter.h": "c",
+        "coefs.h": "c"
+    }
+}

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})

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)
+

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];
+```

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 
+  }
+};

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

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);
+	}
+}
+

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

test/src/IIR_Filter/main.c

@@ -0,0 +1,36 @@
+#include <stdio.h>
+#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;
+}