一、时域自适应线谱增强器

Ruoyo176

已于 2023-11-04 21:26:54 修改

阅读量1.9k

点赞数 7

分类专栏： # 自适应滤波器专栏文章标签： matlab

于 2022-10-09 11:24:59 首次发布

本文链接：https://blog.csdn.net/qq_45490227/article/details/127223438

版权

自适应滤波器专栏专栏收录该内容

3 篇文章 5 订阅

订阅专栏

自适应滤波器

一、时域自适应线谱增强器

文章目录

自适应滤波器
前言
自适应滤波+LMS算法
调用函数
- snr函数：计算信噪比
- SNR_signal函数：计算信噪比
主函数
仿真结果
总结

前言

自适应滤波器是能够根据输入信号自动调整性能进行数字信号处理的数字滤波器。

自适应滤波+LMS算法

function [yn,W,en] = adaptiveFilter_LMS(xn,dn,M,mu,itr)
%自适应滤波器+LMS(Least Mean Squre)算法
%   xn   -->  带噪信号(行向量)
%   dn   -->  原始信号(行向量)
%   M    -->  滤波器阶数(标量)
%   mu   -->  收敛因子(标量)   要求大于0,小于xn的相关矩阵最大特征值的倒数
%   itr  -->  迭代次数(标量)
%   yn   <--  滤波信号(行向量)
%   W    <--  权值矢量(列向量)
%   en   <--  误差(列向量)
% 参数个数必须为4个或5个
if nargin == 4                 % 4个时递归迭代的次数为xn的长度/nargin函数输入参数数目 
    itr = length(xn);          %最大数组维度的长度
elseif nargin == 5             % 5个时满足M<itr<length(xn)
    if itr>length(xn) || itr<M %||或
        error('迭代次数过大或过小!');
    end
else
    error('请检查输入参数的个数!');
end

% 初始化参数
en = zeros(itr,1);             % 误差序列,en(k)表示第k次迭代时预期输出与实际输入的误差/zero创建全零数组
W  = zeros(itr,M);             % 每一行代表一个加权参量,每一列代表-次迭代,初始为0
yn = xn';

% 迭代计算
for k = M:itr                               % 第k次迭代/用来重复指定次数的 for 循环
    x = xn(k:-1:k-M+1);                     % 滤波器M个抽头的输入
    yn(k) = W(k-1,:) * x';                     % 滤波器的输出
    en(k) = dn(k) - yn(k);                % 第k次迭代的误差
    if(mean(abs(en(k-10:k)))<1e-5)
        W(k,:) = W(k-1,:);                  % 滤波器权值计算的迭代式
    else
        W(k,:) = W(k-1,:) + 2*mu*en(k)*x;   % 滤波器权值计算的迭代式
    end
end
yn = yn';
end

调用函数

snr函数：计算信噪比

snr_out = snr(s_k,x_k);   % s_k为参考信号，x_k为输入信号

SNR_signal函数：计算信噪比

自定义的计算信噪比的公式。

function snr = SNR_signal(signal,noiseSignal)
%计算信噪比
%   s_k    --> 原始信号(行向量)
%   x_k    --> 带噪信号(行向量)
%   snr    <-- 信噪比(标)
signalPower = mean(signal.*signal);                          % 求出信号功率
noisePower = mean((noiseSignal-signal).*(noiseSignal-signal));    % 求出噪声功率
snr = 10*log10(signalPower/noisePower);                      % 由信噪比定义求出信噪比，单位为db
end

主函数

%%
close all;clear;clc;
% 1. 创建仿真信号
fs = 16000;             % 采样频率
f0 = [80,160];          % 信号基频
A = [0.5,0.5];            % 幅值
phi = [0,0];            % 初始相位
T = 10;                 % 信号时长
SNR_in = -25;           % 输入信噪比
t = 0:1/fs:T-1/fs;   
s_k = sum(A'.*sin(2*pi*f0'*t+phi'),1);              % 输出原始信号(行向量)       
n_k = randn(size(s_k));                             % 用randn函数产生高斯白噪声
signal_power = mean(s_k.*s_k);                      % 求出信号的平均能量
noise_power = mean(n_k.*n_k);                       % 求出噪声的平均能量
noise_variance = signal_power / ( 10^(SNR_in/10) );    % 计算出噪声设定的方差值
n_k = sqrt(noise_variance/noise_power)*n_k;         % 按噪声的平均能量构成相应的白噪声
x_k = s_k+n_k;                                      % 构成带噪语音

SNR_signal(s_k,x_k)
% 2. 滤波，降采样
s_k_d = lowpass(s_k,1000,fs);
s_k_j = resample(s_k_d,1,5);
x_k_d = lowpass(x_k,1000,fs);
x_k_j = resample(x_k_d,1,5);
snr(s_k_j,x_k_j)
fs_j = fs/5;
plot_s1(s_k_j,x_k_j,fs_j);  % 看一下频谱
%% 3. 线谱增强
M = 2048;               % 滤波器阶数
mu = 0.00001;           % LMS迭代步长
[yn,w,en] = adaptiveFilter_LMS(x_k_j,s_k_j,M,mu);
%%
plot_s2(yn,s_k_j);            % 给出一个误差收敛曲线图
plot_s3(yn,x_k_j,s_k_j,fs_j/2); % 给出一个输入输出信噪比对比图
plot_s4(yn,x_k_j,s_k_j,fs_j); % 给出滤波前后短时傅里叶变化
%%
x_f = 0:fs_j/M:fs_j;
W = fft(w(end,:));
figure;plot(x_f(1:M/2),abs(W(1:M/2)));
%%
function plot_s1(s_k_j,x_k_j,fs)
% 看一下滤波后的输入信号和参考信号频谱
figure;
windowLen = 2^fix(log2(fs));
s1 = s_k_j(1:windowLen);
s2 = x_k_j(1:windowLen);
S1 = abs(fft(s1));
S2 = abs(fft(s2));
x_f = 0:fs/windowLen:fs-fs/windowLen;
subplot(211);plot(x_f(1:windowLen/2),S1(1:windowLen/2));xlabel('频率');ylabel('幅值');title('滤波后参考信号频谱');
subplot(212);plot(x_f(1:windowLen/2),S2(1:windowLen/2));xlabel('频率');ylabel('幅值');title('滤波后输入信号频谱');
end

function plot_s2(yn,dn)
% 给出一个误差收敛曲线图
figure;plot(yn-dn);xlabel('采样点');ylabel('误差');title('训练迭代误差');
end

function plot_s3(yn,xn,dn,fs)
% 每隔N点做一次SNR
len = length(yn);
N = floor(len/fs);
snr_in = zeros(1,N+1);
snr_out = zeros(1,N+1);
snr_gain = zeros(1,N+1);
x_t = 0:1:N;
for i=1:N
    index = (i-1)*fs+1:i*fs;
%     snr_in(i+1) = snr(dn(index),xn(index));
%     snr_out(i+1) = snr(dn(index),yn(index));
%     snr_gain(i+1) = snr_out(i+1)-snr_in(i+1);
    snr_in(i+1) = SNR_signal(dn(index),xn(index));
    snr_out(i+1) = SNR_signal(dn(index),yn(index));
    snr_gain(i+1) = snr_out(i+1)-snr_in(i+1);
end
figure;plot(x_t,snr_in,'o-g');hold on;plot(x_t,snr_out,'o-b');hold on;plot(x_t,snr_gain,'o-r');
xlabel('采样时间/s');ylabel('信噪比/dB');title('信噪比对比图');
legend('输入信噪比','输出信噪比','信噪比增益');
end

function plot_s4(yn,xn,dn,fs)
% 给出滤波前后短时傅里叶变化
figure;
subplot(311);plot_lofar(dn,fs,'参考信号LoFar谱图',0,600);
subplot(312);plot_lofar(xn,fs,'输入信号LoFar谱图',0,600);
subplot(313);plot_lofar(yn,fs,'输出信号LoFar谱图',0,600);
end

function plot_lofar(s,fs,str_title,LF,HF)
%绘制lofar谱图(2048点FFT)
%   s          --> 输入时域信号(行向量)
%   fs         --> 采样频率(标量)
%   str_title  --> 标题
if(nargin==2)
    str_title = 'LoFar谱图';
elseif(nargin==3)
    LF = 0;
    HF = 550;
elseif(nargin==4)
    HF = 550;
elseif(nargin<2 || nargin>5)
    error('请检查输入参数的个数!');
end
windowLen = 2048;
x_i_f = max(1,floor(LF/fs*windowLen)):floor(HF/fs*windowLen);
[S, F, T] = spectrogram(s, hanning(windowLen), windowLen/2, windowLen, fs);
S = mapminmax(S',0,1)';
imagesc(F(x_i_f),T,20*log10(abs(S(x_i_f,:)')));colorbar;colormap("jet");set(gca,'YDir','normal');
ylabel('时间 (s)');xlabel('频率 (Hz)');title(str_title);
end

function snr = SNR_signal(signal,noiseSignal)
%计算信噪比
%   s_k    --> 原始信号(行向量)
%   x_k    --> 带噪信号(行向量)
%   snr    <-- 信噪比(标)
signalPower = mean(signal.*signal);                          % 求出信号功率
noisePower = mean((noiseSignal-signal).*(noiseSignal-signal));    % 求出噪声功率
snr = 10*log10(signalPower/noisePower);                      % 由信噪比定义求出信噪比，单位为db
end