matlab构造线性相位FIR滤波器

前言

用生成的一组音频文件举例

一、构造一组声音

模拟钢琴音乐，采用逐渐衰减振荡的正弦波

FFT的频域展示：

源代码：

function sound_fir
    Fs = 1000; % 采样频率
    freq = [200, 230, 260, 290, 320, 350, 380,410, 440,470]; % 频率数组
    rythm = 0.5; % 持续时间，单位：秒
    gap_duration = 0.2; % 间隔时间，单位：秒
    y = []; % 初始化 y

    note_length = numel(freq); % 使用频率的长度

    % 生成时间向量
    x = linspace(0, rythm, floor(Fs * rythm));
    s = zeros(note_length, numel(x)); 

    % 生成每个频率的波形
    for i = 1:note_length
        s(i, :) = sin(2 * pi * freq(i) * x) .* (1 - x / (rythm)); 
    end

    % 合并音频信号，并添加间隔
    for i = 1:note_length
        y = [y, s(i, :)]; % 添加当前音符
        % 添加间隔
        y = [y, zeros(1, round(Fs * gap_duration))]; % 插入0.2秒间隔
    end

    sound(y, Fs); % 播放合成的音频

    % 做FFT频域变换
    N = length(y); % 采样点数
    frequence = fft(y); % 对时域信号进行FFT变换
    % 幅值修正
    F2 = (frequence / N);
    F1 = F2(1:N/2+1); % 选取前半部分
    F1(2:end-1) = 2 * F1(2:end-1);

    % 画图
    f = Fs * (0:(N/2)) / N;

    % 对其中一段做FFT频域变换
    N0 = length(s(1, :)); % 采样点数
    frequence_single = fft(s(1, :)); % 对第一个音符进行FFT变换
    % 幅值修正
    F20 = (frequence_single / N0);
    F10 = F20(1:N0/2+1); % 选取前半部分
    F10(2:end-1) = 2 * F10(2:end-1);
    
    % 绘图
    f0 = Fs * (0:(N0/2)) / N0;

    subplot(221);
    plot(x, s(2, :)); title('单个输入音频的时域信号'); xlabel('时间 (s)'); ylabel('幅度');%可以设置是哪个频率的音频
    
    subplot(222);
    plot(linspace(0, length(y) / Fs, length(y)), y); title('输入音频的时域信号'); xlabel('时间 (s)'); ylabel('幅度');
    
    subplot(223);
    plot(f, abs(F1)); title('全部频域信号'); xlabel('频率 (Hz)'); ylabel('|F1(f)|'); xlim([0 600]);

    audiowrite('D:\matlab\DSP\sound.wav', y, Fs);

    student_id = [0 0 0 0 0 0 0 0 0 0 0 0 0];
    BP_filter = [];

    % FIR的构造
    for i = 1:length(student_id)
        wc = 2 * pi * freq(student_id(i) + 1) / Fs;
        
        % 频域函数
        N1 = 500; % 窗函数长度(阶数)
        f0 = freq(student_id(i) + 1); % 中心频率 (Hz)
        bw = 20;  % 带宽 (Hz)
        % 创建频率轴
        f1 = linspace(0, Fs, N1);
        % 创建带通滤波器的频率响应
        BP_filter_row = (abs(f1 - f0) <= bw / 2);  % 频率响应
        BP_filter = [BP_filter; BP_filter_row];  % 按行添加到 BP_filter 中
    end
    
    HBP = sum(BP_filter);

    subplot(224);
    for i = 1:length(student_id)
        plot(f1, BP_filter(i, :)); title('构造的理想滤波器——频域'); xlabel('频率 (Hz)'); ylabel('BP_filter'); xlim([0 900]); ylim([0 2]);
        hold on; 
    end   

end

二、采用FIR滤波器做频率筛选

这里采用窗函数构造FIR，将频域上的理想窄带信号，根据IDFT转换为时域的函数，再与汉明窗叠加得到最后的系统函数

  %FIR的构造
    hlp=[];
    fliters=[];
    for i=1:length(student_id)
        %时域函数
        wc=2*pi*freq(student_id(i)+1)/Fs;
        N=500;
        nd=(N+1)/2;%群延时为阶数的一半
        f0 = freq(student_id(i)+1); % 中心频率 (Hz)
        bw = 20;  % 带宽 (Hz)

        wc1=(wc+pi*bw/Fs);%连续转离散，尺度归一化
        wc2=(wc-pi*bw/Fs);
        
        n = linspace(0,N, N);  % 创建n轴，假设滤波器长度为Fs
        hlp_row=(sin(wc1 * (n-nd)) - sin(wc2 * (n - nd))) ./ (pi * (n-nd));   
        hlp_row(isnan(hlp_row)) = wc / pi;  % 处理n=nd时的奇点
        hlp=[hlp;hlp_row];
        %窗函数
        hamming_window = 0.54 - 0.46 * cos(2 * pi * (0:N-1) / (N-1));
        h = hlp_row .* hamming_window; % 应用窗函数
        
        fliters = [fliters;h]; % 存储滤波器
    end