对应matlab的带通butter的C代码实现
时间: 2023-08-08 15:12:37 浏览: 154
以下是一个简单的C语言实现,用于模拟MATLAB中的butter函数设计带通滤波器的数字系数:
```c
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define PI 3.141592654
void butterworth_bandpass(int N, double Wn1, double Wn2, double *b, double *a) {
int i, j;
double alpha, theta;
double *p_real, *p_imag, *p_ang, *p_denom1, *p_denom2;
double *p_num_real, *p_num_imag, *p_denom_real, *p_denom_imag;
double *p_tmp_real, *p_tmp_imag;
p_real = (double *)malloc(sizeof(double) * (N / 2 + 1));
p_imag = (double *)malloc(sizeof(double) * (N / 2 + 1));
p_ang = (double *)malloc(sizeof(double) * (N / 2 + 1));
p_denom1 = (double *)malloc(sizeof(double) * (N / 2 + 1));
p_denom2 = (double *)malloc(sizeof(double) * (N / 2 + 1));
p_num_real = (double *)malloc(sizeof(double) * (N / 2 + 1));
p_num_imag = (double *)malloc(sizeof(double) * (N / 2 + 1));
p_denom_real = (double *)malloc(sizeof(double) * (N / 2 + 1));
p_denom_imag = (double *)malloc(sizeof(double) * (N / 2 + 1));
p_tmp_real = (double *)malloc(sizeof(double) * (N / 2 + 1));
p_tmp_imag = (double *)malloc(sizeof(double) * (N / 2 + 1));
alpha = (N - 1) / 2.0;
for (i = 0; i <= N / 2; i++) {
theta = PI * i / N;
p_real[i] = -2 * alpha * sin(theta);
p_imag[i] = 2 * alpha * cos(theta);
p_ang[i] = atan2(p_imag[i], p_real[i]);
p_denom1[i] = 1 + Wn1 * sin(p_ang[i]);
p_denom2[i] = 1 + Wn2 * sin(p_ang[i]);
}
for (i = 0; i <= N / 2; i++) {
p_num_real[i] = cos(p_ang[i]);
p_num_imag[i] = -sin(p_ang[i]);
p_denom_real[i] = p_denom1[i] * p_denom2[i];
p_denom_imag[i] = 0;
}
for (i = 0; i <= N / 2; i++) {
for (j = 0; j <= N / 2; j++) {
if (i != j) {
p_tmp_real[i] += p_num_real[j] * p_denom_real[(i - j + N / 2) % (N / 2 + 1)] - p_num_imag[j] * p_denom_imag[(i - j + N / 2) % (N / 2 + 1)];
p_tmp_imag[i] += p_num_real[j] * p_denom_imag[(i - j + N / 2) % (N / 2 + 1)] + p_num_imag[j] * p_denom_real[(i - j + N / 2) % (N / 2 + 1)];
}
}
b[i] = p_num_real[i] * p_denom_real[i] - p_num_imag[i] * p_denom_imag[i];
a[i] = p_denom_real[i] * p_denom_real[i] + p_denom_imag[i] * p_denom_imag[i];
}
free(p_real);
free(p_imag);
free(p_ang);
free(p_denom1);
free(p_denom2);
free(p_num_real);
free(p_num_imag);
free(p_denom_real);
free(p_denom_imag);
free(p_tmp_real);
free(p_tmp_imag);
}
int main() {
int N = 4;
double Wn1 = 0.2;
double Wn2 = 0.6;
double *b = (double *)malloc(sizeof(double) * (N / 2 + 1));
double *a = (double *)malloc(sizeof(double) * (N / 2 + 1));
butterworth_bandpass(N, Wn1, Wn2, b, a);
printf("b = {");
for (int i = 0; i <= N / 2; i++) {
printf("%f, ", b[i]);
}
printf("}\n");
printf("a = {");
for (int i = 0; i <= N / 2; i++) {
printf("%f, ", a[i]);
}
printf("}\n");
free(b);
free(a);
return 0;
}
```
请注意,此实现仅考虑了带通滤波器的情况。如果需要实现其他类型的滤波器,需要进行修改。
相关推荐
![doc](https://img-home.csdnimg.cn/images/20210720083327.png)
![docx](https://img-home.csdnimg.cn/images/20210720083331.png)
![docx](https://img-home.csdnimg.cn/images/20210720083331.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)