要移值的matlab函数:
h3 = fir1(16,[0.25 0.50]);
C语言版本
#include <iostream>
#include <cmath>
#define PI acos(-1)
double sincEasy(double *x, int len, int index) {
double temp = PI * x[index];
if (temp == 0) {
return 1.0; // sinc(0) = 1
}
return sin(temp) / temp;
}
double* fir1(int lbflen, double Wn[], double lbf[]) {
double alpha = 0.5 * (lbflen - 1);
double *m = new double[lbflen];
double *R_sin = new double[lbflen];
double *L_sin = new double[lbflen];
for (int i = 0; i < lbflen; i++) {
m[i] = i - alpha;
lbf[i] = 0.0;
}
for (int i = 0; i < 2; i += 2) {
double left = Wn[i];
double right = Wn[i + 1];
for (int j = 0; j < lbflen; j++) {
R_sin[j] = right * m[j];
L_sin[j] = left * m[j];
lbf[j] += right * sincEasy(R_sin, lbflen, j) - left * sincEasy(L_sin, lbflen, j);
}
}
// Apply window function
for (int i = 0; i < lbflen; i++) {
double Win = 0.54 - 0.46 * cos(2.0 * PI * i / (lbflen - 1));
lbf[i] *= Win;
}
// Scaling, if needed
double scale_frequency = 0.5 * (Wn[0] + Wn[1]);
double s = 0.0;
for (int i = 0; i < lbflen; i++) {
s += lbf[i] * cos(PI * m[i] * scale_frequency);
}
for (int i = 0; i < lbflen; i++) {
lbf[i] /= s;
}
delete[] m;
delete[] R_sin;
delete[] L_sin;
return lbf;
}
int main() {
int lbn = 16;
double Wn[2] = { 0.25, 0.50 };
double* lbf = new double[lbn + 1];
lbf = fir1(lbn + 1, Wn, lbf);
// Output coefficients for verification
for (int i = 0; i < lbn + 1; i++) {
std::cout << lbf[i] << " ";
}
std::cout << std::endl;
delete[] lbf;
system("pause");
return 0;
}
结果是完全一样的,OK!