用C++实现计算数字滤波器的频率响应

用GCC下的C++实现设定在绘制频率响应时，将在[0, Fs)范围内使用64个点（即频率分辨率）来计算和显示系统的频率响应，然后传递函数的分子系数[0.008 -0.033 0.05 -0.033 0.008]和分母系数[1 2.37 2.7 1.6 0.41]，这些系数定义了一个有理传递函数，最后设置系统的采样频率为1024，用这些系统参数来分析和展示一个离散时间线性时不变（LTI）系统的频率响应的。使用类似MATLAB的freqz函数计算系统的频率响应如下，以便做出幅频和相频响应曲线。

N=64;
num=[0.008 -0.033 0.05 -0.033 0.008];% 分子系数
den=[1 2.37 2.7 1.6 0.41];% 分母系数
Fs=1024;
freqz(num,den,N,Fs);

这段代码的目的是为了分析和可视化一个由给定分子和分母系数定义的离散时间LTI系统的频率响应。在MATLAB的freqz函数中，它使用了离散傅里叶变换（DFT）的一个子集来计算频率响应。以下是C++代码的实现。

#include <iostream>
#include <vector>
#include <cmath>
#include <complex>

using namespace std;

vector<complex<double>> freqz(const vector<double>& num, 
                                       const vector<double>& den, 
                                       int N, int Fs) {
    int N2 = N / 2 + 1; // 由于对称性，我们只需要一半的频率点
    vector<complex<double>> H(N2);
    vector<complex<double>> W(N); // 复数频率点

    // 生成复数频率点（W = exp(-j*2*pi*k/N)）
    for (int k = 0; k < N; ++k) {
        double wk = -2 * M_PI * k / N;
        W[k] = exp(complex<double>(0, wk));
    }

    // 计算频率响应 H(k) = num(1) + num(2)*W^1 + ... + num(n+1)*W^n / (1 + den(2)*W^1 + ... + den(n+1)*W^n)
    for (int k = 0; k < N2; ++k) {
        complex<double> numerator = 0, denominator = 1;
        for (size_t i = 0; i < num.size(); ++i) {
            numerator += num[i] * pow(W[k], i);
        }
        for (size_t i = 1; i < den.size(); ++i) { // 跳过den[0]，因为它总是1
            denominator += den[i] * pow(W[k], i - 1);
        }
        H[k] = numerator / denominator;
    }

    return H; // 返回频率响应以供进一步处理或绘图
}

int main() {
    int N = 64;
    vector<double> num = {0.008, -0.033, 0.05, -0.033, 0.008};
    vector<double> den = {1, 2.37, 2.7, 1.6, 0.41};
    int Fs = 1024; // 采样频率，但在这里对于freqz计算不是必需的

    vector<complex<double>> H = freqz(num, den, N, Fs);

    // 遍历并打印H的每个复数的实部和虚部
    for (const auto& complexNum : H) {
        cout << "(" << complexNum.real() << ", " << complexNum.imag() << ")" << endl;
    }
    return 0;
}

Fs（采样频率）在freqz函数的计算中不是必需的，因为我们只关心相对频率（从0到π的范围内）。但是，如果您需要绘制与采样频率相关的绝对频率，那么您可以使用Fs来转换频率点。