基于MATLAB的音频信号的抽取,混叠和数字限带滤波(附工程源码、设计报告)

这篇具有很好参考价值的文章主要介绍了基于MATLAB的音频信号的抽取,混叠和数字限带滤波(附工程源码、设计报告)。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

设计要求

设计要求:

1、能够从音频文件读取44100Hz采样频率的x(n),可以通过参数设置读取的起始时间和持续时间;

2、调用MATLAB resample函数对x(n)进行抽取得到y1(m) ;

3、直接对x(n)进行抽取,得到y2(m);

4、先卷积滤波,然后抽取得到y3(m);

5、抽取结合多相滤波得到y4(m);

6、分析和对比不同帧片段的y1(m), y2(m), y3(m)和y4(m)的时域和频域信息;

7、分析和对比运算效率;

8、基于多相滤波结构完成插值和分数倍采样频率变换的实现;

9、滤波效果的基于信噪比定量分析。

以上涉及的抽取因子D均为参数可任意设置,滤波器和具体实现结构需要自行设计。

设计原理

整数因子的抽取

满足关系式xD=(m)=x(mD),如果采样速率为fs,那得到的序列xD(n)的采样率为fs/D,抽取后信号不产生混叠的频率范围变为fm<fs/2D。而抽样序列的频谱如下图

基于MATLAB的音频信号的抽取,混叠和数字限带滤波(附工程源码、设计报告)

图1 抽取后混叠示意图

可以看到抽取序列的频谱是原序列频谱D倍展宽后,按2π的整数倍移位并叠加而成,所以需求中当D=30,采样频率变为1470Hz,小于音频信号的最高频率,必然会发生失真。为保证抽取后不丢失信息或减少丢失信息,应在抽取前先用一个数字低通滤波器进行滤波,对应理想数字滤波器的幅度响应应满足

|HD(ejw)|= D, (|w|<π/D); 0,(π/D≤|w|≤π)

而在实际应用中,我们无法实现理想数字低通滤波器,所以我们需要利用窗函数法,自主设计FIR数字低通滤波器,使其满足对应的幅度响应,达到理想低通滤波器的效果。

整数因子插值

满足xI(m)=x(m/I),(m=0,±I,±2I...);0,(其他),得到频谱X(ejwI)内插后信号频谱被压缩了I倍,需要通过低通滤波器消除内插带来的镜像,滤波器满足

|HI(ejw)| = 1, (|w|<π/I); 0,(π/I≤|w|≤π)

采样率的分数倍转换

即通过D抽取和I插值结合来实现,为了保证不丢失信息,先插值再抽取,考虑到抽样与插值的处理过程均需要使用低通滤波器,所以如图分析可以将两个低通滤波器合二为一,低通截止频率取二者之间的较小值。

基于MATLAB的音频信号的抽取,混叠和数字限带滤波(附工程源码、设计报告)

图2 采样率分数倍转换原理图

多项滤波结构

抽取和插值是多速率信号处理系统中的两个最基本的运算,当对运算速度要求相当高。但在抽取模型中,低通滤波器在减速之前实现的,而插值模型,低通滤波器又位于提速之后,对实时处理不利。所以,为了提高运算速率,抽取利用如图1.3的多项滤波器结构,数字滤波器位于抽取器之后,即滤波是降速后进行的,降低了处理速度的要求,提升实时处理能力。

基于MATLAB的音频信号的抽取,混叠和数字限带滤波(附工程源码、设计报告)

图3 抽取的多项滤波器结构

同理,对于插值的多项滤波器结构使滤波器位于内插器前,使滤波在提速前,降低了对处理速度的要求的同时,也提高了运算精度和降低对字长的要求。

总体模块设计

获取音频数据getdata模块

本模块设计采用matlab内置的audioread函数设置采样点的方式来实现特定起始时间和结束时间的取样,然后使用自定义的get_spectrum画出音频原信号的频谱

function [y1,Fs1] = getdata(path,a,b,n)
Fs=44100;
start_=a*Fs;
end_=start_+b*Fs;
sample=[start_,end_];
[y,Fs1]=audioread(path,sample);
y1=y(:,1);
get_spectrum(y,Fs,'采样44100Hz的原信号频谱',n);
end
基于MATLAB的音频信号的抽取,混叠和数字限带滤波(附工程源码、设计报告)

画频谱函数get_spectrum模块

通过fftshift函数将抽样后的频谱,搬移回零频率附近,因为抽样后的频谱会被搬移到采样频率附近,通过搬移回零频率附近能够重现原始信号的频谱

function get_spectrum(y,Fs,name,n)
N=2^nextpow2(length(y));%得到对应长度,加快fft速度
Yk=fft(y,N);
Ykshift=fftshift(Yk);%频谱搬移
f= (-N/2:N/2-1)*(Fs/N);%频率轴范围搬移到零频率为中心
figure(n);plot(f,abs(Ykshift));
title(name);
xlabel('Hz');
end

使用resample函数抽取resample_D模块

采用matlab内置的resample函数进行在抽样,得到有用信号数据y1,然后通过自己写的画频谱函数模块,画出调用resample后的频谱

function y1= resample_D(y,fs2,fs,n)
disp('resample的时间');
tic
y1=resample(y,fs2,fs);
toc%计算resample抽取时间
get_spectrum(y1,fs2,'resample得到的频谱',n);
end
基于MATLAB的音频信号的抽取,混叠和数字限带滤波(附工程源码、设计报告)

直接抽取direct_D模块

通过直接采用matlab里面切片的用法每隔30个点采一个值得到直接抽取的数据y2,然后使用画频谱函数得到频谱

function y1=direct_D(y,Fs2,D,n)
disp('直接抽取所需时间');
tic
y1=y(1:D:length(y));
toc %计算时间
get_spectrum(y1,Fs2,'直接抽取所得频谱',n);
end
基于MATLAB的音频信号的抽取,混叠和数字限带滤波(附工程源码、设计报告)

FIR滤波器生成FIRfilter_hann模块

采用matlab内置的fir1,设计FIR滤波器,窗函数选用汉宁窗hann,此处函数传参中的归一化的截止频率变为原来的1/D,下采样因子D=30。为了满足设计指标,采用汉宁窗长度为660,然后可以通过通过freqz画出滤波器频响曲线

function b=FIRfilter_hann(N,m)
wc=1/m;%归一化频率变为原来的1/D
b=fir1(N-1,wc,hann(N));
end

直接卷积滤波conv_filter模块

直接通过原音频数据y,与滤波器的单位冲激响应卷积得到限带滤波后的信号,此模块中将y与上个模块的滤波器系数卷积得到数据yf,然后通过30步长直接抽取得到数据y3,然后使用画频谱函数得到频谱

function y1= conv_filter(y,Fs2,b,D,n)
disp('直接卷积filter所需时间');
tic
yf=conv(y,b);
y1=yf(1:D:length(y));
toc
get_spectrum(y1,Fs2,'直接卷积filter所得频谱',n);
end
基于MATLAB的音频信号的抽取,混叠和数字限带滤波(附工程源码、设计报告)

多相滤波(先抽取后滤波)polyphase2模块

本模块通过抽取器结构的I型变换,然后通过恒等变换实现,滤波运算在抽取之后,大大提高了运算速度以及运算精度,减少了运算累积误差,提高了实时处理信号能力,通过此多相滤波后的下采样信号y4的频谱

function bb= polyphase2(y,b,D)
y=reshape(y,1,length(y));%将列向量转换为行向量
n_b=length(b);%滤波器长度
len=floor((length(y)+D-1)/D)+round(length(b)/D)-1;
bb=zeros(1,len);
for i =1:D
    y1=[zeros(1,i-1),y,zeros(1,D-i)];%延时
    y2=y1(1:D:length(y1));%抽取
    hd=b(i:D:n_b);
    y3=conv(y2,hd);%滤波
    if length(y3) >= len%累加
        bb=bb+y3(1:len);
    else
        bb=bb+[y3,zeros(1,len-length(y3))];
    end
end
bb=bb(1:ceil(length(y)/D));bb=reshape(bb,length(bb),1);
end
基于MATLAB的音频信号的抽取,混叠和数字限带滤波(附工程源码、设计报告)

画时域波形get_waveform模块

通过此模块将y1,y2,y3,y4的时域波形全部利用plot函数画出来,传参需要信号的离散点起始点数,以及结束点数,得到四个数据的波形

function get_waveform(y1,y2,y3,y4,start_,end_,n)
a=start_:end_;
figure(n);
subplot(4,1,1);
plot(a,y1(start_:end_));
title('resample的时域波形');
subplot(4,1,2);
plot(a,y2(start_:end_));
title('直接抽取的时域波形');
subplot(4,1,3);
plot(a,y3(start_:end_));
title('直接卷积的时域波形');
subplot(4,1,4);
plot(a,y4(start_:end_));
title('多相滤波的时域波形');
end
基于MATLAB的音频信号的抽取,混叠和数字限带滤波(附工程源码、设计报告)

抽取信噪比计算SNR_cal模块

首先需要找出四种音频信号相互之间的时延,然后进行对齐后再进行信噪比计算,计算时延使用了finddelay函数

得到时延之后首先需要对信号进行对齐,即移动11位,然后计算信噪比,代码如下

function y= SNR_cal(y1,y2,y3,y4)
    y1=y1(1:length(y1)-11);y2=y2(1:length(y2)-11);
    y3=y3(12:length(y3));y4=y4(12:length(y4)); %信号对齐
    s=sum(y1.^2);
    sigma21=sum((y2-y1).^2);
    SNR2=10*log10(s/sigma21);
    sigma31=sum((y3-y1).^2);
    SNR3=10*log10(s/sigma31);
    sigma41=sum((y4-y1).^2);
    SNR4=10*log10(s/sigma41);
    y=[SNR2,SNR3,SNR4];
end

多相滤波插值insert_polyphase2模块

还是采用FIR滤波器设计FIRfilter_hann模块生成插值所需要的消除高频镜像滤波器,此时窗函数长度变为600

设置上采样因子I=30,与前面的抽取结构代码相似,只是把抽取换成了补零。通过理论课所学插值器结构设计insert_polyphase2子函数,利用for循环遍历,实现滤波器系数按相分类,分别运算后累加,实现变换,滤波运算后插值,使得系统运算更快,然后使用plot函数,画出了插值后通过FIR滤波器消除高频镜像后的波形如图所示,然后通过画频谱get_spectrum函数得到,插值消除高频镜像后的频谱

function yf= insert_polyphase2(y,b,I)
 y=reshape(y,1,length(y));N=length(b);
 ttiik=(length(y)+round(N/I)-1)*I+I-1;
 yf=zeros(1,ttiik);
 for i =1:I
      hd=b(i:I:N);
      y1=conv(hd,y);%filter
      y2=zeros(1,length(y1)*I);%插值补零
      y2(1:I:length(y2))=y1;
      y2=[zeros(1,i-1),y2,zeros(1,I-i)];%延时
      if length(y2) >= ttiik%累加
        yf=yf+y2(1:ttiik);
     else
        yf=yf+[y2,zeros(1,ttiik-length(y2))];
     end
 end
 yf=yf(1:I*length(y));
 yf=reshape(yf,length(yf),1);% turn it into column vector
end
基于MATLAB的音频信号的抽取,混叠和数字限带滤波(附工程源码、设计报告)

主函数调用联合调试

clc;
clear all;
close all;
%% 实现四种抽取、两种限带滤波、以及时域频域对比分析、信噪比分析
start_=1; end_=5; %参数设置读取的起始时间和持续时间为1-5s
[y,Fs]=getdata('audio.mp3',start_,end_,1); %1.1调用自定义函数,读取音频数据,并画出原频谱图,由start和end1控制读取时间
D=30; %抽取因子D为30
Fs2=Fs/D;
firN_d=660; %有限长滤波器截取长度660
y1=resample_D(y,Fs2,Fs,2); %调用自定义函数resample_D,实现matlab自带下采样限带滤波,并计时
y2=direct_D(y,Fs2,D,3); %1.3调用direct_D函数,直接抽取,并计时
b=FIRfilter_hann(firN_d,D); %创建hanning窗有限长滤波器      
y3 = conv_filter(y,Fs2,b,D,4); %1.4先卷积滤波,然后抽取
disp('多相滤波2所需时间');
tic
y4=polyphase2(y,b,D); %1.5先多相滤波,再抽取
toc
get_spectrum(y4,Fs2,'多相滤波频谱',5);
SNR=SNR_cal(y1,y2,y3,y4); 
get_waveform(y1,y2,y3,y4,1,length(y4),6);
%% 插值的实现_based on polyphase filter algorithm
I=30;
firN_i=600;   
Fs3=Fs*I;
bi=FIRfilter_hann(firN_i,I);
y6=insert_polyphase2(y,bi,I);
n=1:length(y6);
get_spectrum(y6,Fs3,'多相滤波insert频谱',7);axis([-2000 2000 0 4500]);
figure(8);
plot(n,y6);title('多相滤波插值后的波形');
%% 分数倍的速率转换__based on polyphase filter algorithm
I1=35;D1=55;
firN_i1=600;firN_d1=600;
Fs4=floor(Fs*I1/D1);
bi1=FIRfilter_hann(firN_i1,I1);
bd=FIRfilter_hann(firN_d1,D1);
yi=insert_polyphase2(y,bi1,I1);
yd=polyphase2(yi,bd,D1);
figure(9);
n1=1:length(yd);
plot(n1,yd);
title('多相滤波速率变换waveform');
get_spectrum(yd,Fs4,'多相滤波速率变换的频谱',10);

时域分析

四种信号的时域波形,通过肉眼可以观察出,通过了限带滤波的resample、直接卷积滤波、多相滤波的方式中时域波形相差不大,而没有通过限带滤波的直接抽取方式,时域波形明显失真较大;然后通过sound函数播放了进行30抽取后的音频发现直接抽取方式的音频有明显的混叠现象,混叠失真相对其他三种方式非常明显。

频域分析

观察得,原始信号能量频率绝大部分集中在5kHz以内,根据奈奎斯特抽样定理,信号采样率至少大于10kHz,频谱不会发生明显混叠,当信号采样率为44100/30=1470时频谱必然产生混叠现象,抽样后不能还原原始信号,音频必然会有失真。对比四种抽取后频谱可以发现,直接抽取的频谱在700Hz,400Hz附近能够观察出明显的频谱混叠现象

虽然混叠失真无法避免,但是通过前置限带滤波器,将信号频带限制在抽取后的折叠频率内,能够一定程度上减少混叠失真,通过四个频谱可以发现,通过自己设计FIR滤波器起到了较为良好的抗混叠效果

运算效率分析

对比发现直接抽取的运算时间最短,这是符合逻辑的,因为直接抽取只是进行了简单的按步长进行抽取,而没有通过滤波运算,所以所需运算时间是最短的。然后发现直接卷积滤波的时间要远远小于多相滤波的时间,这与理论课中所学知识是不符的,因为我这里采用了多相滤波2,是滤波器I型变换后的等效结构,滤波运算是在抽取降速之后进行的,按理论来说,运算速度应该快于直接卷积滤波,这里的原因,我猜测是因为我在进行子函数polyphase2设计时,采用了for循环进行按相分类,执行for循环时,matlab只能挨着按顺序进行遍历,最后把运算结果相加,浪费掉了大量运算资源,在实际中,滤波器的多相结构各支路运算应是同时进行的,最后全部加起来。所以由于for循环的特性,导致了多相滤波运算效率大大减小,运算时间远远大于直接卷积。

信噪比指标分析

调用信噪比计算模块计算SNR2,SNR3,SNR4,三个信噪比结果分别为-0.7185dB,29.84dB,29.84dB,然后计算得到SNR3-SNR2=30.5585dB>30dB,|SNR4-SNR3|=0<1dB,满足设计指标。

通过分析可以看到,直接抽取得到的信号信噪比非常低甚至小于0,说明直接抽取信号失真较大,而通过了限带滤波器后再抽取的信号增大了约30dB,所以限带滤波器大幅提高了信号信噪比,在一定程度上减小了混叠的影响

分数倍速率变换分析

有理论课中学到,为了避免丢失信息,应先进行插值后进行抽取,这里设置下采样因子D=55,上采样因子I=35,然后结合前面所做的polyphase2,和insert_polyphase2模块很容易得出速率变换为以前的I/D=7/11后的时域波形与频谱

基于MATLAB的音频信号的抽取,混叠和数字限带滤波(附工程源码、设计报告)

此时采样率变为floor(44100*7/11) = 28063Hz通过分析对比以及结合理论知识,可以看到,时域采样率提高的同时,在频域上频谱也被进行了压缩,同时我调用sound函数播放了速率变换后的音频,发现与原音频相差并不大,这是因为原音频频谱能量绝大部分都集中于5kHz内,而变换后采样频率28063Hz远远大于2fm,根据奈奎斯特抽样定理,音频没有明显的混叠失真文章来源地址https://www.toymoban.com/news/detail-411335.html

如有任何疑问可以留言

声明:只用作分享学习,共同进步,其他不做任何用途

到了这里,关于基于MATLAB的音频信号的抽取,混叠和数字限带滤波(附工程源码、设计报告)的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包