精华内容
下载资源
问答
  • MATLAB滤波程序

    2018-07-03 20:54:31
    一段信号滤波程序MATLAB实现的,不仅仅一种滤波啊
  • 心电信号处理matlab程序对学习、课程设计以及毕业设计有比较大的帮助
  • Matlab语音信号滤波程序
  • 心电信号matlab滤波

    2018-04-16 10:21:43
    程序是利用matlab设计切比雪夫滤波器实现心电信号滤波,并绘制滤波前后的时域波形。
  • Matlab语音信号滤波程序.
  • matlab进行信号滤波

    2014-08-15 11:13:00
    matlab进行信号滤波,附带程序,几个特别经常用的例子,希望大家有帮助
  • 基于matlab的语音信号滤波处理

    千次阅读 2021-08-25 11:54:02
    基于matlab的语音信号滤波处理摘要:本课程设计的主要目的是在MATLAB环境下,使用窗口设计法设计一个滤波器,并语音信号进行滤波去噪。开发平台为MATLAB,设计方法为窗口设计法。用麦克风采集一段语音信号,绘制...

    基于matlab的语音信号滤波处理

    摘要:本课程设计的主要目的是在MATLAB环境下,使用窗口设计法设计一个滤波器,并对语音信号进行滤波去噪。开发平台为MATLAB,设计方法为窗口设计法。用麦克风采集一段语音信号,绘制波形并观察其频谱,给定相应技术指标,用凯塞窗设计一个满足指标的FIR滤波器,对该语音信号进行滤波去噪处理,比较滤波前后的波形和频谱并进行分析。在滤波前后语音信号波形的变化中,由于我们录制的语音信号噪声不大,所以观察并不明显,但在频域波形中,我们可以明显的看到设计的滤波器对语音信号进行了滤波处理,将噪声进行了滤除。此次滤波基本达到了要求,完成了设计指标。


    关键词:滤波去噪;FIR滤波器;凯塞窗;MATLAB

    目录


    1 课程设计研究内容 1

    1.1课程设计研究内容 1

    1.2课程设计步骤及流程图 1

    1.3课程设计要求 3

    2 设计原理 3

    2.1 MATLAB简介 3

    2.2 滤波器 3

    3 设计与实现过程 3

    3.1实现方法 3

    3.1.1 FIR滤波器 4

    3.1.2窗口设计法 4

    3.1.3凯塞窗 4

    3.2 设计过程 4

    3.2.1录制语音信号 4

    3.2.2 对原信号加入噪声 5

    3.2.2 滤波器的设计 7

    3.2.2信号的滤波处理 9

    3.3 仿真结果与分析 10

    4 调试问题与解决方案 11

    5总结与展望 11

    致谢 12

    参考文献 12

    附录 12

    附录1. 滤波器脉冲响应源程序 12

    附录2. 理想低通滤波器计算源程序 13


    1 课程设计研究内容

    1.1课程设计研究内容

    1.语音信号的采集

    在Windows下录制一段格式为.wav的语音,利用函数wavread对语音信号进行采样。

    2.语音信号的频谱分析

    用MATLAB程序对原始语音信号进行采样、频谱分析,并绘制出采样后语言信号时域波形图和频谱图,并针对此图分析语音信号特点。

    3.语音信号加噪与频谱分析

    利用MATLAB程序产生信号噪声,并加入到语音信号中,模仿语音信号被污染,并对其频谱分析,与原始语音信号进行对比,分析差异。

    4.设计数字滤波器

    根据语音信号的特点,设计数字滤波器,对加噪后的语音信号进行滤波处理。

    5.验证滤波器的滤波效果

    对滤波后的语音信号进行时域、频域分析,并将滤波前后的时域波形、频谱波形进行相比较,分析信号的变化,从而验证所设计滤波器的滤波效果是否达到了滤除高频噪音、保留低频原始语音信号的目的。

    6.回放语音信号

    利用函数sound对滤波后语音信号进行回放。

    1.2课程设计步骤及流程图

    设计本课题的流程为:采集一段语音信号。将语音信号的文件名命名为input18.wav,再用MATLAB中的wavread函数求出语音信号的三个参数,分别为:每个样本的值,生成该语音波形文件时的采样频率,波形文件样本的码数,再对信号及加入单频干扰后的语音信号做傅立叶变化,绘制出时域和频域的波形。最后通过滤波绘制滤波前后时域波形对比图和幅频特性对比图,并回放滤波前后的语音信号来验证是否达到去噪的目的。课程的设计流程图如图1-1所示:

    v2-7a9b092650ad672b7451651b86931c86_b.jpg
    图1-1设计流程图


    2 设计原理

    2.1 MATLAB简介

    MATLAB是一个为科学和工程计算专门设计的交互式大型软件,是一个可以完成各种精确计算和数据处理的,可视化的,强大的计算工具。它具有丰富的函数资源和工具箱资源,语言精练,代码灵活,面向对象,控制功能优良,图形工能也强大。并且它的兼容性很好,几乎能在所有的PC机和大型计算机上运行,适用于Windows,UNIX和多种系统平台。MATLAB形形色色的工具箱中包括控制系统,信号处理,小波分析,统计,优化等,能够很好的运用于语音信号的滤波去噪。

    2.2 滤波器

    数字滤波器在数字信号处理的各种应用中发挥着十分重要的作用。它是通过对采样数据信号进行数学运算处理来达到滤波的目的。FIR数字滤波器设计的方法有三类,一类是窗口设计法(时间窗口法),第二类是频率采样法,第三类是等波纹优化设计。时间窗口设计法是从单位脉冲响应序列着手,使h(n)逼近理想的单位脉冲响应序列hd(n);频率采样法是使所设计的FIR数字滤波器的频率特性某些离散频率点上的值准确地等于所需滤波器在这些频率点处的值,在其它频率处的特性则有较好的逼近。等波纹优化设计也叫最佳一致逼近准则,最佳一致逼近即选择N个频率采样值(或时域h(n) 值),在给定频带范围内使频响的最大逼近误差达到最小。可保证局部频率点的性能也是最优的,误差分布均匀,相同指标下,可用最少的阶数达到最佳化。本次课程设计采用的就是窗口设计法。


    3 设计与实现过程

    3.1实现方法

    3.1.1 FIR滤波器

    数字滤波器(Digital Filter,简称为DF)是指用来对输入信号进行滤波的硬件和软件。 所谓数字滤波器,是指输入、输出均为数字信号,通过一定运算关系改变输入信号所含频率成分的器件。数字滤波器和模拟滤波器相比,因为信号的形式和实现滤波的方法不同,数字滤波器具有比模拟滤波器精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配等优点。一般用两种方法来实现数字滤波器:一是采用通用计算机,把滤波器所要完成的运算编成程序通过计算机来执行,也就是采用计算机软件来实现;二是采用实际专用的数字处理硬件。

    数字滤波器按照冲激响应的时域特性可分为:无限长单位冲响应滤波器(IIR)和有限长单位冲击响应滤波器(FIR),但与IIR相比,在满足同样阻带衰减的情况下需要的阶数较高,滤波器的阶数越高,占用的运算时间就越多,因此在满足指标要求的情况下应尽量减少滤波器的阶数。

    FIR滤波器的基本结构可以理解为一个分节的延时线,把每一节的输出加权累加,可得到滤波器的输出,FIR滤波器的冲激响应h(n)是有限长的,数学M阶FIR滤波器可以表示为:

    v2-1a096c97cc6abd28ace60188589dc491_b.jpg


    (3-1)


    v2-b8b6475df7e246f398ebc461848e94bb_b.jpg

    (3-2)


    3.1.2窗口设计法

    窗口设计法的基本思想是要选取某一种合适的理想频率选择性滤波器(这种滤波器总是有一个非因果,无限长的脉冲响应),然后将它的脉冲响应截断(或加窗)以得到一个线性相位和因果的FIR滤波器。

    3.1.3凯塞窗

    窗函数的主瓣宽度和旁瓣峰值衰耗是矛盾的,一项指标的提高总是以另一项指标的下降为代价,窗口选择实际上是对两项指标作权衡。而两项指标是跳变的,于是有人提出可调整窗,适当修改参数,可在这两项指标间作连续的选择。常用的可调整窗是凯塞(Kaiser)窗。凯塞(Kaiser)窗全面地反映主瓣与旁瓣衰减之间的交换关系, 可以在它们两者之间自由地选择它们的比重。

    凯塞窗的表达式是:

    v2-24d26d60dc39af44c75ec5ba4e9e5217_b.jpg


    v2-98d5fff4c206173e156faf402bc3bfb9_b.png


    (3-3)


    式中,

    v2-b2366c37962184a72d86178a260b6fb8_b.png

    是第一类修正的零阶贝塞尔函数,可以用级数展开来计算它的值。

    3.2 设计过程

    3.2.1录制语音信号

    用windows工具中的录音机录制一段语音信号,语音为“请鞭挞我吧,公瑾!”,时间长度约为2s。将语音信号的文件名设置为input18.wav,并将文件保存在MATLAB下的WOK文件夹里面。然后在MATLAB平台上,用wavread函数调出此语音信号,并得到其采样率fs和比特数bits。

    [x,fs,bits]=wavread('input18.wav') ; % 输入参数为文件的全路径和文件名,输出的第一个参数是信号的样本值,fs是生成该波形文件时的采样率,bits是波形文件每样本的编码位数。

    sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放录音

    x=x(:,2);

    N=length(x); % 计算信号x的长度

    fn=1000; % 单频噪声频率

    t=0:1/fs:(N-1)/fs; % 计算时间范围,样本数除以采样频率

    x=x';

    y=x+0.1.*sin(fn*2*pi*t); %y为加入单频干扰信号后的语音

    sound(y,fs,bits) ;

    在MATLAB平台上,用plot函数画出原始语音信号,如图3-1所示:


    v2-0dccbb9dcefe8c2e50900c0dc73f3f03_b.jpg



    图3-1


    3.2.2 对原信号加入噪声

    在MATLAB平台上,对原始信号和加噪信号进行fft变换,取幅度谱,并对频谱进行分析。具体实现如下:

    [x,fs,bits]=wavread('input18.wav') ; % 输入参数为文件的全路径和文件名,输出的第一个参数是信号的样本值,fs是生成该波形文件时的采样率,bits是波形文件每样本的编码位数。

    sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放录音

    x=x(:,2);

    N=length(x); % 计算信号x的长度

    fn=1000; % 单频噪声频率

    t=0:1/fs:(N-1)/fs; % 计算时间范围,样本数除以采样频率

    x=x';

    y=x+0.1.*sin(fn*2*pi*t); %y为加入单频干扰信号后的语音

    sound(y,fs,bits) ; % 对加噪信号进行回放

    X=abs(fft(x));% 对原始信号进行fft变换,取幅度谱

    Y=abs(fft(y)); % 对加噪后信号进行fft变换,取幅度谱

    X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分

    deltaf=fs/N; % 计算频谱的谱线间隔

    f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围


    subplot(2,2,1);plot(t,x); axis([0 1 -1 1]); grid on; %画原始语音信号的时域图

    xlabel('时间(s)');ylabel('幅度');title('原始语音信号');

    subplot(2,2,2);plot(f,X);axis([0 4000 0 1500]); grid on; %画原始语音信号幅度谱图

    xlabel('频率(Hz)');ylabel('幅度谱');title('语音信号幅度谱图');

    subplot(2,2,3);plot(t,y); axis([0 1 -1 1]); grid on; %画加噪语音信号的时域图

    xlabel('时间(s)');ylabel('幅度');title('加入单频干扰后的语音信号');

    subplot(2,2,4);plot(f,Y);axis([0 4000 0 1500]); grid on; %画加噪语音信号幅度谱图

    xlabel('频率(Hz)');ylabel('幅度谱');title('加入干扰后的语音信号幅度谱图');

    在MATLAB上运行结果如图3-2所示:


    v2-56281a2ebd67efd9994873f3f46314be_b.jpg


    图3-2原始信号及加噪信号的时域图及幅度谱图

    由图3-2可以看出,对比原始信号与加噪信号幅度谱图我们可以清晰的看到在1000Hz是出现了一个脉冲,从而使原有的音频听上去有了噪声,这样就到达了预期结果

    3.2.2 滤波器的设计

    滤波器设计就是要找到一组能满足特定滤波要求的系数向量a和b,而它主要是通过设计指标来实现的。滤波器设计的要求或指标一般是在频域上给出的,常用的滤波器频域指标有:fp1、fs1、fs2、fp2、Rp、As。要达到最佳的滤波效果,则需要对fp1、fs1、fs2、fp2和Rp、As进行适当的调整。由图3-3可以看出,语音信号可以选择fp1=900;fs1=950;fs2=1050;fp2=1100;Rp=1;As=60的滤波器。在MATLAB中,通常采用1/2采样频率进行归一化处理,如果将频率转化为角频率,则需将归一化频率乘以pi。设计程序如下:

    fp1=900;fs1=950;fs2=1050;fp2=1100;Rp=1;As=60; % 带阻滤波器设计指标

    fc1=(fp1+fs1)/2;fc2=(fp2+fs2)/2; df=min((fs1-fp1),(fp2-fs2)); % 计算上下边带中心频率,和频率间隔

    wc1=fc1/fs*2*pi;wc2=fc2/fs*2*pi; dw=df/fs*2*pi; %将Hz为单位的模拟频率换算为rad为单位的数字频率

    ws1=fs1/fs*2*pi;ws2=fs2/fs*2*pi;

    M=ceil((As-7.95)/(2.285* dw)+1)+1; % 计算汉宁窗设计该滤波器时需要的阶数

    n=0:M-1; beta=0.1102*(As-8.7); % 定义时间范围

    w_kaiser= kaiser(M,beta); % 产生M阶的凯塞窗

    hd_bs=ideal_lp(wc1,M)+ideal_lp(pi,M)-ideal_lp(wc2,M); % 调用自编函数计算理想带阻滤波器的脉冲响应

    h_bs=w_kaiser'.*hd_bs; % 用窗口法计算实际滤波器脉冲响应

    [db,mag,pha,grd,w]=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性

    Rp=-min(db(wc1/dw+1:1:wc2/dw)); As=max(-round(db(ws2/dw+1:1:501)));

    subplot(2,2,1);plot(w/pi,db);axis([0 0.2 -70 10]);

    xlabel('w/pi');ylabel('dB');title('滤波器幅度响应图');grid on;

    subplot(2,2,2);plot(w/pi,mag);axis([0 0.2 0 1.2]);

    xlabel('w/pi');ylabel('幅度mag');title('滤波器幅度响应图');grid on;

    subplot(2,2,3);plot(w/pi,pha);axis([0 1 -3 3]);

    xlabel('w/pi');ylabel('相位pha');title('滤波器相位响应图');grid on;

    subplot(2,2,4);stem(n,h_bs);axis([0 1500 0 1]);

    xlabel('n');ylabel('h(n)');title('滤波器脉冲响应图');grid on;

    所设计的滤波器的幅度响应和脉冲响应如图3-3所示:


    v2-9d1bac662fbbb3b9d4a906d1ffe76c97_b.jpg


    图3-3滤波器的频率特性图



    由3-3图可以看出Rp接近于0,而As大于开始的设定值(As=60db),由此可见满足了设计要求,达到了设计目的。

    3.2.2信号的滤波处理

    滤波器设计完成后,在MATLAB平台上调用函数fftfilt对信号实现滤波,绘制语音信号去噪前后时域图,频域幅度图形并进行比较。具体程序如下:

    y_fil=fftfilt (h_bs, y); % 用设计好的滤波器对y进行滤波

    Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半

    subplot(3,2,1);plot(t,x); axis([0 3 -0.2 0.2]);

    xlabel('时间(s)');ylabel('幅度');title('原始语音信号x');grid on;

    subplot(3,2,2);plot(f,X);axis([0 4000 0 70]);

    xlabel('频率(Hz)');ylabel('幅度谱');title('语音信号幅度谱X');grid on;

    subplot(3,2,3);plot(t,y); axis([0 3 -0.2 0.2]);

    xlabel('时间(s)');ylabel('幅度');title('加干扰后的语音信号x1');grid on;

    subplot(3,2,4);plot(f,Y);axis([0 4000 0 70]);

    xlabel('频率(Hz)');ylabel('幅度谱');title('加干扰语音信号幅度谱X1');grid on;

    subplot(3,2,5);plot(t,y_fil);axis([0 3 -0.2 0.2])

    xlabel('时间(s)');ylabel('幅度');title('滤波后语音信号y');grid on;

    subplot(3,2,6);plot(f,Y_fil); axis([0 4000 0 70])

    xlabel('频率(Hz)');ylabel('幅度谱');title('滤波后语音信号幅度谱Y');grid on;

    sound (y_fil,fs,bits)


    在MATLAB平台上的运行结果如图3-4所示:


    v2-7935d6e507d0e0aca4d238b5ace32d5a_b.jpg


    图3-4-语音信号去噪前后时域图及频域幅度图


    由图3-4可以看出,在原始信号上加入的1000Hz 的脉冲信号已经被设计的滤波器滤去,符合实验所要到达的预期效果。

    3.3 仿真结果与分析

    通过观察滤波前后语音信号波形的变化,即观察图3-3,我们可以知道,原信号的时域图与滤波去噪信号的时域图基本相似;原信号与滤波去噪信号的频谱图波形也大致相似。通过观察可以看到,加噪信号的时域图中大部分都被加入的噪声给遮盖了,加噪信号的频谱图中,我们可以很明显地看到与原信号频谱图相比,它在频率1000Hz的地 方有一个尖脉冲,而滤波去噪信号的频谱图中该尖脉冲已经消失,幅度变为了零,波形大致与原图相似,可见滤波去噪效果基本不错。在将三个信号的时域波形和频谱图比较之后,我们还要通过回放去滤波去噪音乐信号,来跟原信号相比,以检验滤波器的效果。在MATLAB中,函数sound可以对声音进行回放。其调用格式为:sound (x,fs,bits)。用sound(y_fil,fs,bits)语句回放该滤波去噪信号,便可以感觉到滤波后的语音号与原信号差不多。


    4 调试问题与解决方案

    1.在采用凯塞窗口设计的FIR滤波器时得不到理想的滤波器;

    解决方法:这是因为通带截止频率与阻带起始频率之间的差值设置的太小或者太大,通过适当的选择参数fp1、fs1、fs2及fp2,绘制出来的图形效果比较明显,基本符合设计指标。这样通过MATLAB运算出来的滤波器的阻带波纹达到了滤波器设计的要求,得到了比较理想的滤波器。

    2.设计滤波器绘制出来的频谱图中不知道As和Rp是否达标,根本不了解在什么样的情况下As和Rp才能符合滤波器的要求;

    解决方法:通过请教老师和寻找同学的帮助解决这一问题,慢慢的得到改善,最终符合设计要求。

    3. 在滤波的过程中对噪声的幅度值滤波滤不干净,在噪声处总有一些幅度值;

    解决方法:反复观察程序的错误,进行修改,和同学交流,最终得到改善。



    附录

    附录1. 滤波器脉冲响应源程序

    function [db,mag,pha,grd,w] = freqz_m(b,a); 
    % freqz 子程序的改进版本
     
    % ------------------------------------ 
    % [db,mag,pha,grd,w] = freqz_m(b,a);  
    % db = [0 到pi弧度]区间内的相对振幅(db) 
    % mag = [0 到pi弧度]区间内的绝对振幅
    % pha = [0 到pi弧度]区间内的相位响应
    % grd = [0 到pi弧度]区间内的群迟延
    % w = [0 到pi弧度]区间内501个频率样本向量
    % b = Ha(z)的分子多项式系数(对FIR b=h) 
    % a = Ha(z)的分母多项式系数(对FIR: a=[1]) 
    %  
    [H,w] = freqz(b,a,1000,'whole');  
    H = (H(1:1:501))'; w = (w(1:1:501))'; 
    mag = abs(H);  
    db = 20*log10((mag+eps)/max(mag)); 
    pha = angle(H); 
    % pha = unwrap(angle(H)); 
    grd = grpdelay(b,a,w); 
    % grd = diff(pha); 
    % grd = [grd(1) grd];  
    % grd = [0 grd(1:1:500); grd; grd(2:1:501) 0]; 
    % grd = median(grd)*500/pi; 



    附录2. 理想低通滤波器计算源程序

    function hd = ideal_lp(wc,M); 
    % 理想低通滤波器计算
     
    % -------------------------------- 
    % [hd] = ideal_lp(wc,M)  
    % hd = 0 to M-1之间的理想脉冲响应
    % wc = 截止频率(弧度) 
    % M = 理想滤波器的长度
    % 
    alpha = (M-1)/2; 
    n = [0:1:(M-1)]; 
    m = n - alpha + eps; 
    hd = sin(wc*m) ./ (pi*m);

    展开全文
  • MATLAB处理语音信号滤波程序.doc
  • matlab低通滤波程序

    2013-07-30 10:16:37
    %实现了频率为20和200Hz单频叠加cos信号的低通滤波,使输出仅含有20Hz分量 clear; fs=1200; %采样频率 N=1200; % N/fs 秒数据 n=0:N-1; t=n/fs; %时间 fL=20; fH=200; sL=cos(2*pi*fL*t); sH=cos(2*pi*fH*t); ...
  • 使用MATLAB进行卡尔曼滤波(9轴)IMU数据 这是用于9轴IMU传感器的卡尔曼滤波算法。 (加速度计,陀螺仪,磁力计) 您可以看到带有数据的图形化IMU传感器。 测验 示范 -将很快添加。 特征 动画情节 时间线 硬铁偏置...
  • 维纳滤波matlab代码语音信号“实时声音平台” 您可以在这个平台上做的事情:•看声音•听声音•测量声音•修改声音 该平台是可以以任何明智的方式连接的模块的集合。 例如:录制声音,其进行操作,其进行可视化...
  • 对于心电信号处理不错的资源,本次实验是要我们先产生心电信号,再加噪声再滤波。由一个心电信号数据表输出心电信号,由于已知的心电信号有噪声,所以我先将信号进行滤波,得到正确的心电信号。再加噪声,通过滤波...
  • 表面肌电信号处理的matlab程序,包括带通滤波、50Hz陷波滤波程序,以及计算时域、频域的指标iMEG、RMS , MF、MPF
  • Matlab自带函数的滤波程序,各种滤波方法都可用
  • 应用Matlab对人体的心电信号进行滤波实验目的综合应用信号频谱分析和数字滤波器设计的知识,实现心电信号滤波。加深理解信号时域和频域分析的物理概念,理解设计指标的工程概念,认识不同类型滤波器的特性和适用...

    应用Matlab对人体的心电信号进行滤波

    实验目的

    综合应用信号频谱分析和数字滤波器设计的知识,实现心电信号的滤波。加深理解信号时域和频域分析的物理概念,理解设计指标的工程概念,认识不同类型滤波器的特性和适用范围。

    实验环境

    1.微型电子计算机(PC);

    2.安装Windows10操作系统,MATLAB等开发工具。

    实验原理

    首先对待滤波的心电信号进行频谱分析,观察信号频率分布的规律,从而确定数字滤波器的类型(FIR滤波器、IIR滤波器、自适应滤波器、小波滤波器等)。在加性噪声的情况下,若信号的频谱与噪声的频谱基本不重叠,可以采用频率选择滤波器(FIR滤波器、IIR滤波器)。若信号的频谱与噪声的频谱重叠较多,可以采用自适应滤波、小波滤波等。若为乘性噪声,可以根据同态滤波的原理对信号进行预处理,然后再按照加性噪声的情况处理。

    在确定了数字滤波器的类型后,还需要根据信号时域特性、频域特性、或时频特性确定滤波器的设计参数,设计出相应的数字滤波器。

    最后,利用该数字滤波器对信号进行滤波,在时域和频域观察信号滤波的主观及客观效果。若主观及客观效果满足要求,说明分析过程和滤波方法正确有效,若不满足要求,需要重新分析和设计。

    实验内容和任务要求

    人体的心电信号通常分布在200Hz的范围内,在测量过程中往往会受到工业高频噪声的干扰,所以必须经过低通滤波处理后,才能作为判断心脏功能的有用信息。若已知一个实际心电信号的采样序列样本如下:

    x(n)={-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0}

    其中存在高频干扰。要求:

    (1)设计一个合适的滤波器,对上述心电信号的采样序列进行滤波处理,滤除其中的

    干扰成分,画出滤波器的幅频响应和相频响应曲线。

    (2)分别显示滤波前后心电信号的时域波形和幅度谱,观察总结滤波作用与效果。

    展开全文
  • 基于MATLAB的语音滤波实验实验目的:在Matlab环境下语音的频谱进行处理(数字滤波)并试听效果;在Matlab环境下语音的抽样率进行处理(语音压缩)并试听效果实验步骤:一、音频文件的压缩(抽取)。利用windows...

    基于MATLAB的语音滤波实验

    实验目的:

    1. 在Matlab环境下对语音的频谱进行处理(数字滤波)并试听效果;
    2. 在Matlab环境下对语音的抽样率进行处理(语音压缩)并试听效果

    实验步骤:

    一、音频文件的压缩(抽取)。

    1. 利用windows附件中的录音机功能录制8~10秒的.wav语音文件,并以lei为文件名保存到Matlab/work的文件夹中。

    a.打开 开始/程序/附件/娱乐/录音机;


    v2-33bf7d0199153a332d4411c2be34b462_b.jpg


    b.用windows media player播放一首音乐并用MIC对着耳机录音或自已说话录音(按

    v2-329e9b841d608da8fe53a8f4e0e0de1d_b.jpg

    键),到10秒时停止(按

    v2-975742155984f0381c9d1cba3c3e77b7_b.jpg

    键);

    c.将录制的文件加存为C:/Matlab/work中,文件名为leii.wav;


    v2-b827015fcd6efad16c6e0b93d3a35255_b.jpg

    v2-08f44185aa6a631ea0848932bdc4de28_b.jpg

    v2-f0ae34482b30a9ca5c0b75d4e20eef5f_b.jpg



    1. 打开Matlab并新建一.m文件;
    2. 在.m文件中用y=wavread(‘lei.wav’)命令读入语音文件。


    v2-d03e1b30e637cd3b363a5c3ef586feab_b.jpg


    1. 语音压缩:在m命令窗中输入如下命令:

    v2-f8a1b7aed11353f13412e8a9ddca2a43_b.jpg


    1. 运行sample2.m之后会在work文件夹中生成一个名为lei2的.wav文件,如下图:


    v2-e25f83b4f9ce403c1e5bb118cc8def9e_b.jpg


    1. 双击lei2音频文件,用耳机试听效果,并跟lei1的效果比较。
    2. 在sample2.m文件中改变抽取倍率s (必须为正整数),重复4、5、6步,观察在不同抽取倍率s下的音频质量,

    (注意:在运行sample2.m之前必须将work中名为lei2的.wav音频文件删除,或在.m文件中wavwrite()中的保存文件名改为其它的名字。)

    二、音频信号的时域滤波(音频数据的时域卷积)。

    (一)、低通滤波

    1. 打开Matlab并新建一.m文件,在.m文件中用y=wavread(‘lei.wav’)命令读入语音文件。
    2. 在m命令窗中输入如下命令,并加存为sample3.m,运行该m文件。


    v2-2f001b56c8f88706a21795fb6c6a9ba6_b.jpg


    1. 双击lei3音频文件,用耳机试听效果,并跟lei1的效果比较。
    2. 再加一级h(n)的低通滤波,重复2、3步,如下图:

    (注意:在运行lei2.m之前必须将work中名为lei3的.wav音频文件删除,或在.m文件中wavwrite()中的保存文件名改为其它的名字。)


    v2-f6e0127948e98f1baec7fd099f3465b2_b.jpg


    1. 重复2、3、4步,观察在不同阶数的低通滤波下的音频质量。

    (二)、高通滤波

    1.打开Matlab并新建一.m文件,在.m文件中用y=wavread(‘lei.wav’)命令读入语音文件。

    2.在m命令窗中输入如下命令,并加存为sample4.m,运行该m文件。


    v2-6744a88a056a36fde37492f68d6031e0_b.jpg


    1. 双击lei4音频文件,用耳机试听效果,并跟lei1的效果比较。
    2. 再加一级h(n)的低通滤波,重复2、3步,如下图:

    (注意:在运行lei2.m之前必须将work中名为lei4的.wav音频文件删除,或在.m文件中wavwrite()中的保存文件名改为其它的名字。)


    v2-33804f180842bacd95ef1877d2f9a0bc_b.jpg


    再加一级高通滤波:


    v2-d84de39573a5ea0fd43ea345932a0774_b.jpg


    (三)时域低通滤波时频域的频谱变化:

    1.打开在第(一)步中创建的sample3.m文件,并在原文件中加入以下命令,另存为sample33.m


    v2-8a15aace5a0c4df781506eed14cdcdec_b.jpg


    频谱如下图:


    v2-572f9213d3eb85b84ede2cd3f889e392_b.jpg


    2.下图为h(n)为31点的三重低通滤波程序和频谱图:


    v2-3cf7673e7635764f25ff17d33a242126_b.jpg



    v2-9a08a761ade1805da38e126c1b1d401a_b.jpg




    (四)时域高通滤波时频域的频谱变化:

    1.打开在第(二)步中创建的sample4.m文件,并在原文件中加入以下命令,另存为sample44.m


    v2-0ddde8d979c41910d893837b953b9d09_b.jpg


    频谱如下图:


    v2-975e7344e2c29c2c30016ca6f970d1db_b.jpg



    2.下图为h(n)为3点的三重低通滤波程序和频谱图:


    v2-54e3c0f32c1266571016ef603290a713_b.jpg



    v2-d01a352737631f81afde5e3a55785536_b.jpg



    二、音频信号的频域理想滤波处理:

    原音频信号的抽样频率:

    v2-34ddc2b22db9a43ae5b25f0951e79f20_b.jpg

    ,


    v2-cd0e95d0055f09d9d44fa28fc05e0ab5_b.jpg



    v2-0389918b96c5765b2ec1ef14b084b916_b.jpg

    ,该式即为模拟域频率f跟频率域(FFT变换)中k的关系

    (一)理想低通滤波:

    1.设计一截止频率为

    v2-680b9b875e6cfb03a699ea3f5f33f038_b.jpg


    对应的数字频率域(FFT)的系统函数的频率特性为:

    v2-01cb4596f184be251efeb3ca34281768_b.jpg


    2.按下图所示输入m文件,另存为lowfilter1.m ,并运行该程序。


    v2-d0dd078ccb1b578fc6e36360ecf585ae_b.jpg


    2.双击lei5音频文件,用耳机试听效果,并跟lei1的效果比较。

    3.将该m文件中的f0=2000分别改为1000、500、300、200、4000、…….后,运行程序试听效果。

    4.在该m 文件中加入如下命令,重复第3步,并观察Xw、Hw、Yw的频谱。


    v2-a96405a602b0f360c580f99ecd70cc96_b.jpg



    v2-6667cf4652ef26972d0f0537159d57e3_b.jpg





    (二)理想高通滤波:

    1.设计一截止频率为

    v2-9cd36ab899c52f74b7daf23cd85f91d8_b.jpg


    对应的数字频率域(FFT)的系统函数的频率特性为:

    v2-12158b818c300ee7006c8dd50c4fe7a2_b.jpg


    2.按下图所示输入m文件,另存为lowfilter1.m ,并运行该程序。


    v2-4efb32c068e2e77cf38d613c82867d5f_b.jpg


    2.双击lei6音频文件,用耳机试听效果,并跟lei1的效果比较。

    3.将该m文件中的f0=1000分别改为800、500、300、200、1500、2000…….后,运行程序试听效果。

    4.在该m 文件中加入如下命令,重复第3步,并观察Xw、Hw、Yw的频谱。


    v2-69785e78b869c200ff592eac5ec15660_b.jpg



    v2-1f192dfb0a52413f121d18f858a206bb_b.jpg


    (二)理想高通滤波:

    1.设计一截止频率为

    v2-43191b4308b5a20a784b814e3fdab14a_b.jpg


    对应的数字频率域(FFT)的系统函数的频率特性为:

    v2-5d6c25852d30f50abe3557c84c787697_b.jpg


    2.按下图所示输入m文件,另存为lowfilter1.m ,并运行该程序。


    v2-7c37e55c6631d841be295b40be431b10_b.jpg


    2.双击lei7音频文件,用耳机试听效果,并跟lei1的效果比较。

    3.将该m文件中的fl和fh分别改.后,运行程序试听效果。

    4.在该m 文件中加入如下命令,重复第3步,并观察Xw、Hw、Yw的频谱。


    v2-e382b1b6ae4fc9af53756791d1e60dc7_b.jpg

    v2-80e6b35b130449a1ef1ff4659eb63d6a_b.jpg


    基于MATLAB的语音加去噪和延时混响实验

    实验说明:

    1.本实验提供的beiguo.wav,lei1.wav,music.wav,shao.wav,wang.wav均为原始语音信号.

    2.本实验中的jiazao.m为语音加噪实验,xiaozao.m为语音消噪实验,musicadd.m为语音全成实验,musicfilter.m为语音滤波实验,dlaymusic.m为语音混响实验.

    实验步骤:

    1.将本文件夹中的所有.m文件和所有原始语音信号都复制到MATLAB的work文件夹中。

    2.打开MATLAB程序。

    一、语音消噪实验


    %%%%%%在语音中加噪声%%%%%%%

    x1=wavread('lei1.wav');%读取原语音信号,lei1中无噪声.

    fs=22050; %原语音信的采样率为22050Hz

    fn=1000; %设定噪声的频率为1000Hz

    t=1:length(x1); %设置噪声的度度跟原语音信一样长,

    x2=2*sin(2*pi*fn/fs*t);%产生幅度为2频率为fn的正弦波作为噪声.

    x=x1+x2'; %将原子核语音信号跟噪声相加,x为带有噪声的语音信号.

    wavwrite(x,22050,'lei2.wav');%将带有噪声的语音信号转换为声音,lei2中将有噪声


    1).在MATLAB中打开名为jiazao.m的程序,运行该程序,将在work中产生一个新的语音文件lei2.wav

    2).通过试听对比lei1.wav和lei2.wav语音,看噪声是否加上。


    二、语音消噪实验


    %%%%%%消除语音中的噪声%%%%%%%

    x1=wavread('lei2.wav');%读取原语音信号,lei2中带有噪声.

    y=filter(hn,1,x); %将带有噪声的语音信号x经过带阻滤波器进行滤波,以达到消噪目的.

    %Bndstop,FIR,Equiripple,Minimum order,Fs=22050,Fpass1=950,Fstop1=980,Fstop2=1020,Fpass2=1050,Apass1=1,Astop=60,Apass2=1

    wavwrite(y,22050,'lei3.wav');%将经带阻滤波消噪后的信号转换为语音,lei3中将不再有噪声


    1).在MATLAB中打开名为xiaozao.m的程序(暂时不运行)

    2).在MATLAB左下角start中打开FDATool界面,按本程序m文件中注释的参数设计带阻滤波器,

    并通过File-Export-填hn,将设计的滤波器系数导到工作空间。

    3).运行该程序,将在work中产生另一个新的语音文件lei3.wav

    4).通过试听对比lei2.wav和lei3.wav语音,看噪声是否消除。


    三、语音滤波实验


    %%%%%带阻和低通滤波%%%%%%%

    x=wavread('shao.wav');%读取名为shao.wav的原语音信号

    y=filter(hns,1,x); %带阻滤波,滤波器在FDATool中设计,并导到工空间,因本人的中低音太重,高音不足

    %Hpass,FIR,Equiripple,Minmum order,Fs=22050,Fpass1=100,Fstop1=1500,Fstop2=1600,Fpass2=3000,Apass1=1,Astop=30,Apass=1.

    yy=filter(hnh,1,y); %高通滤波,滤波器在FDATool中设计,并导到工空间,因本人的低音太重,高音不足

    %Hpass,FIR,Equiripple,Minmum order,Fs=22050,Fstop=10,Fpass=4000,Astop=20,Apass=1

    wavwrite(yy,22050,'shao2.wav');%将经混响后的信转换为语音,shao2.wav的语音中的中低频分量将有所衰减.


    1).在MATLAB中打开名为musicfilter.m的程序(暂时不运行)

    2)按本程序m文件中注释的参数分别设计带阻和高通滤波器,并通过File-Export-填hns和hnh,将设计的滤波器系数导到工作空间。

    3).运行该程序,将在work中产生另一个新的语音文件lei4.wav

    4).通过试听对比shao.wav和shao2.wav语音,看语音有可不同。


    四、语音混响实验


    %%%%%延时混响%%%%%%%

    x=wavread('shao2.wav'); %读入原始声音

    n=1200; %设定延迟时间t=n/fs秒,改变该数据可改变混响深度(时间间隔)

    N=60; %y设定延迟级数为N级,改变该数据可改变次数

    x1=[x;zeros(N*n,1)]; %将x通过补零延长到经N级延时后的长度

    for i=1:N %进行N次延时,第一次延时在x前补n 个0,后补(N-1)*n个0

    x2=[zeros(i*n,1);x;zeros((N-i)*n,1)]; %第i次延时在x前补i*n个0,后补(N-i)*n个0

    x1=x1+1/(2*i)*x2; %将经延时的信号x1跟x逐次相加

    end

    wavwrite(x1,22050,'shao3.wav'); %将混响后的数据转换为声音


    1).在MATLAB中打开名为dlaymusic.m的程序

    2).运行该程序,将在work中产生另一个新的语音文件shao3.wav

    3).通过试听对比shao.wav、shao2.wav和shao3.wav语音,看语音有可不同。

    4).修改本程序中的n和N,并重复2)和3)的步骤。


    五、语音合成实验


    %%%%%将两首语音全成一首%%%%%%

    m1=wavread('beiguo.wav');%读取一首语音m1

    m2=wavread('wang.wav'); %读取另一首语音m2

    if length(m1)>length(m2) %比较两首语音的长度,将短的补成跟长的相等

    m3=[m2;zeros((length(m1)-length(m2)),1)];

    else

    m3=[m1;zeros((length(m2)-length(m1)),1)];

    end

    m=0.6*m1+m3; %将两个语音相加,为分辩明,将其中一个衰减

    wavwrite(m,22050,'mu3.wav'); %将合成后的信转为语音


    1).在MATLAB中打开名为musicadd.m的程序

    2).运行该程序,将在work中产生另一个新的语音文件mu3.wav

    3).通过试听对比beiguo.wav、wang.wav和mu3.wav语音,看语音有否不同。

    展开全文
  • matlab仿真心电图滤波程序

    热门讨论 2009-11-12 22:48:45
    本例是一个应用savitzky-Golay滤波器滤除心电图噪音的应用例子。
  • matlab文件以LFM信号为例详细介绍了信号匹配滤波的仿真方法和实现,包括时域方法和频域方法
  • 一维中值滤波函数,调用函数,输入滤波窗口值和初始一维滤波向量,输出滤波后的向量
  • 基于matlab心电信号ECG滤波处理 二、源代码 function varargout = ECG(varargin) % ECG MATLAB code for ECG.fig % ECG, by itself, creates a new ECG or raises the existing % singleton*. % % H = ECG returns ...

    一、心电信号数字滤波处理简介

    心电信号作为一种人体的基本生理信号, 是心脏电活动在人体体表的表现, 信号一般比较微弱, 频率在0.05Hz~100Hz范围内, 幅度为10V (胎儿) ~5m V (成人) , 心电信号信噪比和频率都较低, 在心电的采集、放大、检测等过程中, 易受到外界的各种干扰。常见的噪声干扰有:第一是基线漂移, 一般是由人体呼吸和心肌兴奋所引起的, 它的频率低于0.5Hz, 属于低频干扰;其次是肌电干扰, 它是由人体肌肉颤动所致, 它的发生频率具有随机性, 范围在5Hz~2000Hz之间;第三是工频干扰, 它是由室内照明及动力设备影响到人体的分布电容所引起的, 频率为50Hz。消除或减少这些干扰时识别心电信号特征和参数的前提。心电信号噪声来源不同, 频率也存在差异, 正是由于这些差异, 对不同的信号干扰其滤波方法也不同, 滤波可以用硬件实现, 但实现过程相对困难, 也可以用软件编程方法实现, 数字滤波技术成为目前滤除心电干扰的有效手段。

    1 程序设计与实现
    “心电信号的数字滤波处理”软件所要实现的功能和任务如下所示:
    (1) 信号输入:信号源的读取及参数的输入;
    (2) 信号滤波:选择信号分析通道, 选择滤波器类型和种类, 根据选择的滤波器类型及参数指标, 用相对应的阶数选择函数返回阶数N及截止频率Wn, 根据N及Wn利用IIR滤波器响应的设计函数对信号进行滤波处理。
    (3) 滤波器特性演示:显示所设计的各个滤波器的幅度和相位响应;
    (4) 信号显示:包括原始信号的显示和经过各次滤波后信号的显示。
    在这里插入图片描述
    图1 程序设计流程图

    2 具体界面设计如下所示
    (1) 打开并选择文件:创建打开文件对话框并显示文件存储路径, 数据文件通常为.txt或.dat格式。“选择文件”用按钮 (push button) 实现, 当点击时, 能够打开如图4所示的对话框, 可选择数据文件, 并在文本框 (edit text) 内可显示文件存储路径。

    (2) 信号通道选择:由于所采集的心电信号数据是12通道的, 进行心电信号分析时只需选择其中之一, 信号选择通道用下拉菜单 (pop-up menu) 实现。

    (3) 选择滤波器类型:滤波器类型共有四种:Butterworth、Chebyshev1、Chebyshev2、Elliptic, 用下拉菜单 (pop-up menu) 实现。

    (4) 滤波器功能实现:带阻、高通、低通分别用三个单选按钮 (radio button) 实现, 并用按钮组 (button group) 把三个控件组织在同一区域内。

    (5) 参数输入:采样频率 (Fs) 、通带截止频率 (Fp1, Fp2) 、阻带截止频率 (Fs1, Fs2) 、通带波动 (Rp) 、阻带衰减 (Rs) 从界面上输入, 显示这些参数指标的组件为文本标签 (static text) , 显示输入参数的组件为文本框 (edit text) , 当选择带阻时能显示全部参数指标, 当选择低通、高通时, 通过设置属性“visible”, 隐藏通带截止频率 (Fp2) 、阻带截止频率 (Fs2) , 使其编辑框不可见。
    (6) 滤波结果显示:滤波前后的波形在相应的坐标轴上显示, 其中第一个坐标轴显示滤波前的原始心电信号signal, 第二个坐标轴显示signal经过带阻滤波后的信号波形s1, 第三个坐标轴显示s1经高通滤波后的信号波形s2, 第四个坐标轴显示s2经过低通滤波后的波形s3, 第五个坐标轴同时显示原始信号signal和经过三次滤波后的波形s3, 以便进行对比分析, 点击每个坐标轴旁边的“Run”按钮可显示坐标轴相对应的信号波形。

    (7) 滤波器性能演示及退出界面:用按钮 (push button) 实现, 当点击“滤波器性能演示”时, 能够显示选择的滤波器的幅度及相位响应, 当点击“Quit”时, 退出当前窗口, 返回编辑界面。

    二、部分源代码

    function varargout = ECG(varargin)
    % ECG MATLAB code for ECG.fig
    %      ECG, by itself, creates a new ECG or raises the existing
    %      singleton*.
    %
    %      H = ECG returns the handle to a new ECG or the handle to
    %      the existing singleton*.
    %
    %      ECG('CALLBACK',hObject,eventData,handles,...) calls the local
    %      function named CALLBACK in ECG.M with the given input arguments.
    %
    %      ECG('Property','Value',...) creates a new ECG or raises the
    %      existing singleton*.  Starting from the left, property value pairs are
    %      applied to the GUI before ECG_OpeningFcn gets called.  An
    %      unrecognized property name or invalid value makes property application
    %      stop.  All inputs are passed to ECG_OpeningFcn via varargin.
    %
    %      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
    %      instance to run (singleton)".
    %
    % See also: GUIDE, GUIDATA, GUIHANDLES
    
    % Edit the above text to modify the response to help ECG
    
    % Last Modified by GUIDE v2.5 22-May-2020 21:35:12
    
    % Begin initialization code - DO NOT EDIT
    gui_Singleton = 1;
    gui_State = struct('gui_Name',       mfilename, ...
                       'gui_Singleton',  gui_Singleton, ...
                       'gui_OpeningFcn', @ECG_OpeningFcn, ...
                       'gui_OutputFcn',  @ECG_OutputFcn, ...
                       'gui_LayoutFcn',  [] , ...
                       'gui_Callback',   []);
    if nargin && ischar(varargin{1})
        gui_State.gui_Callback = str2func(varargin{1});
    end
    
    if nargout
        [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
    else
        gui_mainfcn(gui_State, varargin{:});
    end
    % End initialization code - DO NOT EDIT
    
    
    % --- Executes just before ECG is made visible.
    function ECG_OpeningFcn(hObject, eventdata, handles, varargin)
    % This function has no output args, see OutputFcn.
    % hObject    handle to figure
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    % varargin   command line arguments to ECG (see VARARGIN)
    
    % Choose default command line output for ECG
    handles.output = hObject;
    
    % Update handles structure
    guidata(hObject, handles);
    
    % UIWAIT makes ECG wait for user response (see UIRESUME)
    % uiwait(handles.figure1);
    
    
    % --- Outputs from this function are returned to the command line.
    function varargout = ECG_OutputFcn(hObject, eventdata, handles) 
    % varargout  cell array for returning output args (see VARARGOUT);
    % hObject    handle to figure
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    
    % Get default command line output from handles structure
    varargout{1} = handles.output;
    
    
    
    function edit1_Callback(hObject, eventdata, handles)
    % hObject    handle to edit1 (see GCBO)
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    
    % Hints: get(hObject,'String') returns contents of edit1 as text
    %        str2double(get(hObject,'String')) returns contents of edit1 as a double
    
    
    % --- Executes during object creation, after setting all properties.
    function edit1_CreateFcn(hObject, eventdata, handles)
    % hObject    handle to edit1 (see GCBO)
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    empty - handles not created until after all CreateFcns called
    
    % Hint: edit controls usually have a white background on Windows.
    %       See ISPC and COMPUTER.
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
    
    
    % --- Executes on button press in pushbutton1.
    function pushbutton1_Callback(hObject, eventdata, handles)
    % hObject    handle to pushbutton1 (see GCBO)
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    [FileName PathName]=uigetfile('*.mat','MAT-files(*.mat)','选择文件');
    str=[PathName FileName];
    set(handles.edit1,'string',str);
    Fs=200;N=512;MEAN=0;MIN=0;MAX=0;VAR=0;STD=0;RR=0;
    global im;global MEAN;global MIN;global MAX;global VAR;global STD;global RR;
    if strcmp(str,'D:\心电数据\被试2  心电\chenwei1.mat')==1
        im=1;
    else
        im=2;
    end
    qq=get(handles.popupmenu1,'value')
    bb=get(handles.popupmenu2,'value')
    switch(im)
        case 1
                load('chenwei1.mat');
                ECG=xin(qq,:);
        case 2
               load('fanglipeng1.mat');
                ECG=xin(qq,:);
    end
    
    % uu1=get(handles.radiobutton7,'value');
    % uu2=get(handles.radiobutton8,'value');
    % uu3=get(handles.radiobutton8,'value');
    % if uu1==1
    %     MAX=max(EGC);
    %     set(handles.edit3,'string',MAX);
    % end
    % if uu1==2
    %     MIN=min(EGC);
    %     set(handles.edit4,'string',MIN);
    % end
    % if uu1==3
    %     MEAN=mean(EGC); 
    %     set(handles.edit5,'string',MEAN);
    % end
    if          bb==1
                disp([num2str(bb)]);
                t=1/50:1/50:length(ECG)/50;
                MEAN=mean(xin(qq,:));MAX=max(xin(qq,:)); MIN=min(xin(qq,:));
                VAR=var(xin(qq,:),0,2);STD=std(xin(qq,:),0,2);
                axes(handles.axes3);
                grid on;plot(t,ECG);
                xlabel('时间(s)');ylabel('幅值(V)');
                title('原始心电信号时域显示');
                N=512;
                y=fft(ECG,N);
                mag=abs(y);
                f=(0:length(y)-1)*Fs/length(y);
                axes(handles.axes5);
                plot(f,mag)
                xlim([0,100]);
                title('原始心电信号频谱图')
                xlabel('频率(Hz)');ylabel('幅值');
                axes(handles.axes8);
                histogram(ECG);xlabel('信号幅值');ylabel('个数');
                title('不同区间心电信号的分布');
    elseif      bb==2
                disp([num2str(bb)]);
                %-----------------------------心电信号高通滤波器进行基线漂移纠正--------------------------------%
                %还可采用椭圆滤波器或中值法(myMedfilt)等进行基线漂移纠正。
                fp=1.4;fr=0.05;rp=1;rs=30;
                wp=fp/(Fs/2);wr=fr/(Fs/2);
                [n,wn]=buttord(wp,wr,rp,rs);
                [b,a]=butter(n,wn,'high');
                [hw,w]=freqz(b,a);
                y2=filter(b,a,ECG);
                axes(handles.axes1);
                grid on;
                t1=0:1/50:(length(y2)-1)/50;
                plot(t1,y2)
                title('基线漂移去除时域')
                xlabel('时间(s)');ylabel('幅值(V)');
                y2=fft(y2,N);
                mfb=abs(y2);
                f=(0:length(y2)-1)*Fs/length(y2);
                axes(handles.axes2);
                grid on;
                plot(f,mfb);
                xlim([0,100]);
                title('基线纠正后信号幅频')
                xlabel('频率(Hz)');ylabel('幅值');
    elseif      bb==3
                disp([num2str(bb)]);
                fp=60;fs=100;                    %通带截止频率, 阻带截止频率  
                rp=1.4;rs=1.6;                   %通带、阻带衰减  
                wp=2*pi*fp;ws=2*pi*fs;     
                [n,wn]=buttord(wp,ws,rp,rs,'s'); %’s’是确定巴特沃斯模拟滤波器阶次和3dB  截止模拟频率  
                [z,P,k]=buttap(n);   %设计归一化巴特沃斯模拟低通滤波器,z为极点,p为零点和k为增益  
                [bp,ap]=zp2tf(z,P,k);  %转换为Ha(p),bp为分子系数,ap为分母系数  
                [bs,as]=lp2lp(bp,ap,wp); %Ha(p)转换为低通Ha(s)并去归一化,bs为分子系数,as为分母系数  
                [hs,ws]=freqs(bs,as);         %模拟滤波器的幅频响应  
                [bz,az]=bilinear(bs,as,Fs);     %对模拟滤波器双线性变换  
                [h1,w1]=freqz(bz,az);         %数字滤波器的幅频响应  
                m=filter(bz,az,ECG);
                t1=0:1/50:(length(m)-1)/50;
                axes(handles.axes1);
                grid on;axis tight;
                plot(t1,m)
                title('低通滤波心电信号时域图');
                xlabel('时间(s)');ylabel('幅值(V)');
                y1=fft(m,N);
                mfa=abs(y1);
                f=(0:length(y1)-1)*Fs/length(y1);
                axes(handles.axes2);
                grid on;axis tight;
                plot(f,mfa);
                xlim([0,100]);
                title('低通滤波心电信号频谱图');
                xlabel('频率(Hz)');ylabel('幅值');
    elseif     bb==4 
               disp([num2str(bb)]);
               %50Hz陷波器:由一个低通滤波器加上一个高通滤波器组成  
                %而高通滤波器由一个全通滤波器减去一个低通滤波器构成  
                Me=100;               %滤波器阶数  
                L=100;                %窗口长度  
                beta=100;             %衰减系数  
                Fs=200;  
                wc1=49/Fs*pi;     %wc1为高通滤波器截止频率,对应51Hz  
                wc2=50/Fs*pi     ;%wc2为低通滤波器截止频率,对应49Hz  
                h=ideal_lp(0.132*pi,Me)-ideal_lp(wc1,Me)+ideal_lp(wc2,Me); %h为陷波器冲击响应  
                w=kaiser(L,beta);  
                y=h.*rot90(w);         %y为50Hz陷波器冲击响应序列  
                y3=filter(y,1,ECG); 
                t1=0:1/50:(length(y3)-1)/50;
                axes(handles.axes1);
                grid on;axis tight;
                plot(t1,y3)
                title('去除工频干扰时域时域')
                xlabel('时间(s)');ylabel('幅值(V)');
                y3=fft(y3,N);
                mfc=abs(y3);
                f=(0:length(y3)-1)*Fs/length(y3);
                axes(handles.axes2);
                grid on;axis tight;
                plot(f,mfc);
                xlim([0,100]);
                title('去除工频干扰时域频域');
                xlabel('频率(Hz)');ylabel('幅值');
    end
    %%
    %心率计算
    %低通滤波
    fp=60;fs=100;                    %通带截止频率,阻带截止频率  
    rp=1.4;rs=1.6;                   %通带、阻带衰减  
    wp=2*pi*fp;ws=2*pi*fs;     
    [n,wn]=buttord(wp,ws,rp,rs,'s'); %’s’是确定巴特沃斯模拟滤波器阶次和3dB  截止模拟频率  
    [z,P,k]=buttap(n);   %设计归一化巴特沃斯模拟低通滤波器,z为极点,p为零点和k为增益  
    [bp,ap]=zp2tf(z,P,k);  %转换为Ha(p),bp为分子系数,ap为分母系数  
    [bs,as]=lp2lp(bp,ap,wp); %Ha(p)转换为低通Ha(s)并去归一化,bs为分子系数,as为分母系数  
    [hs,ws]=freqs(bs,as);         %模拟滤波器的幅频响应  
    [bz,az]=bilinear(bs,as,Fs);     %对模拟滤波器双线性变换 
    [h1,w1]=freqz(bz,az);         %数字滤波器的幅频响应  
    m=filter(bz,az,ECG);
    fp=1.4;fr=0.05;rp=1;rs=30;
    wp=fp/(fs/2);wr=fr/(fs/2);
    [n,wn]=buttord(wp,wr,rp,rs);
    [b,a]=butter(n,wn,'high');
    [hw,w]=freqz(b,a);
    y2=filter(b,a,m);
    %高通滤波
    fp=1.4;fr=0.05;rp=1;rs=30;
    wp=fp/(Fs/2);wr=fr/(Fs/2);
    [n,wn]=buttord(wp,wr,rp,rs);
    [b,a]=butter(n,wn,'high');
    [hw,w]=freqz(b,a);
    y2=filter(b,a,m);
    %带阻滤波
    Me=100;               %滤波器阶数  
    L=100;                %窗口长度  
    beta=100;             %衰减系数  
    wc1=49/Fs*pi;     %wc1为高通滤波器截止频率,对应51Hz  
    wc2=50/Fs*pi     ;%wc2为低通滤波器截止频率,对应49Hz  
    h=ideal_lp(0.132*pi,Me)-ideal_lp(wc1,Me)+ideal_lp(wc2,Me); %h为陷波器冲击响应  
    w=kaiser(L,beta);  
    y=h.*rot90(w);         %y为50Hz陷波器冲击响应序列  
    y3=filter(y,1,y2); 
    %------------------------R波检测及心率计算---------------------------------
    a=1;b=[-0.000563925539225382,-0.000687497480673608,-0.000400105459896818,1.40107308696153e-18,0.000167430218323421,-4.18540736387490e-19,-0.000215512204999462,-2.31378858508372e-19,0.000825087119402599,0.00173867121147225,0.00170041274003132,-2.89185229056478e-18,-0.00281247259984511,-0.00484192899664537,-0.00404159503274541,6.02125250728051e-18,0.00516947529617596,0.00792307431385663,0.00588976752279114,-7.41129393462523e-18,-0.00583817525454588,-0.00764316928687033,-0.00463943332803296,-9.61260744539564e-19,0.00196521978920439,-2.40802770032636e-18,-0.00240551785915901,4.71696771800100e-18,0.00853358817557591,0.0172851206658474,0.0162833499451329,-1.24036605494920e-17,-0.0253107164699268,-0.0426277900469897,-0.0350794664744733,1.73330580157676e-17,0.0449057473944898,0.0701489638765231,0.0540262508439817,-3.98392673063995e-17,-0.0616816640708770,-0.0915638898806427,-0.0672043113722621,-2.75356818775531e-17,0.0700978687033445,0.0996538686301192,0.0700978687033445,-2.75356818775531e-17,-0.0672043113722621,-0.0915638898806427,-0.0616816640708770,-3.98392673063995e-17,0.0540262508439817,0.0701489638765231,0.0449057473944898,1.73330580157676e-17,-0.0350794664744733,-0.0426277900469897,-0.0253107164699268,-1.24036605494920e-17,0.0162833499451329,0.0172851206658474,0.00853358817557591,4.71696771800100e-18,-0.00240551785915901,-2.40802770032636e-18,0.00196521978920439,-9.61260744539564e-19,-0.00463943332803296,-0.00764316928687033,-0.00583817525454588,-7.41129393462523e-18,0.00588976752279114,0.00792307431385663,0.00516947529617596,6.02125250728051e-18,-0.00404159503274541,-0.00484192899664537,-0.00281247259984511,-2.89185229056478e-18,0.00170041274003132,0.00173867121147225,0.000825087119402599,-2.31378858508372e-19,-0.000215512204999462,-4.18540736387490e-19,0.000167430218323421,1.40107308696153e-18,-0.000400105459896818,-0.000687497480673608,-0.000563925539225382];
    d=filter(b,a,y3);
    t=(0:(length(d)-1))/Fs;
    d1=sort(d,'descend');d2=0;
    for i=201 :2200
        d2=d1(i)+d2;
    end;
    d3=d2/2000;
    [R_V,R_L]=findpeaks(d(12001:24000),'minpeakdistance',1,'minpeakheight',d3);
    d3=d2/2000;
    [R_V,R_L]=findpeaks(d(12001:24000),'minpeakdistance',1,'minpeakheight',d3);
    %根据位置和采样频率来计算采样区间内的平均心率
    H_Rate = 60*(length(R_L)-1)/((R_L(length(R_L))-R_L(1))/200);
    RR=num2str(H_Rate,3);
    %算出采样的时间区间
    %R_Time = R_L(length(R_L))/200; 
    
    
    % --- Executes on button press in radiobutton7.
    
    function radiobutton7_Callback(hObject, eventdata, handles)
    % hObject    handle to radiobutton7 (see GCBO)
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    global MAX;
     set(handles.edit3,'string',MAX);
    
    % Hint: get(hObject,'Value') returns toggle state of radiobutton7
    
    
    % --- Executes on button press in radiobutton8.
    function radiobutton8_Callback(hObject, eventdata, handles)
    % hObject    handle to radiobutton8 (see GCBO)
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    global MIN;
    set(handles.edit4,'string',MIN);
    % Hint: get(hObject,'Value') returns toggle state of radiobutton8
    
    
    % --- Executes on button press in radiobutton9.
    function radiobutton9_Callback(hObject, eventdata, handles)
    % hObject    handle to radiobutton9 (see GCBO)
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    global MEAN;
     set(handles.edit5,'string',MEAN);
    % Hint: get(hObject,'Value') returns toggle state of radiobutton9
    
    
    
    function edit3_Callback(hObject, eventdata, handles)
    % hObject    handle to edit3 (see GCBO)
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    
    % Hints: get(hObject,'String') returns contents of edit3 as text
    %        str2double(get(hObject,'String')) returns contents of edit3 as a double
    
    
    % --- Executes during object creation, after setting all properties.
    function edit3_CreateFcn(hObject, eventdata, handles)
    % hObject    handle to edit3 (see GCBO)
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    empty - handles not created until after all CreateFcns called
    
    % Hint: edit controls usually have a white background on Windows.
    %       See ISPC and COMPUTER.
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
    
    

    三、运行结果

    在这里插入图片描述

    四、matlab版本及参考文献

    1 matlab版本
    2014a

    2 参考文献
    [1] 沈再阳.精通MATLAB信号处理[M].清华大学出版社,2015.
    [2]高宝建,彭进业,王琳,潘建寿.信号与系统——使用MATLAB分析与实现[M].清华大学出版社,2020.
    [3]王文光,魏少明,任欣.信号处理与系统分析的MATLAB实现[M].电子工业出版社,2018.
    [4]董兵,超于毅,李振新.基于MATLAB的心电信号的数字滤波处理[J].数字技术与应用. 2012,(10)

    展开全文
  • 表面肌电信号处理的matlab程序,包括带通滤波、50Hz陷波滤波程序,以及计算时域、频域的指标iMEG、RMS , MF、MPF
  • matlab 七种滤波程序

    2014-03-13 21:24:11
    该文件含有matlab编写的7中信号滤波程序
  • 多相滤波MATLAB程序

    热门讨论 2014-05-07 21:28:39
    MATLAB多相滤波的仿真程序,包含子通道的中心频率的设计,解释详细,很快就能看懂。
  • 详细的心电信号噪声处理,包含各种滤波器及小波变换,傅里叶变换等的心电信号处理程序代码。提供完整心电信号数据,UI界面操作,下载即可使用。
  • Matlab 实验报告(题目二)(题目二)声音信号的采集与滤波处理(采用IIR滤波器或FIR滤波器)参考资料:信号的采集、数字信号处理及滤波实例要求:(1)采集声音信号或打开已录好的声音文件,并显示其信号图与频域图。...
  • 一维二维中值滤波均值滤波matlab编码实现
  • 输入信号,加入高斯白噪声,使用维纳滤波进行消除高斯白噪声,得要期望的信号,显示输入。噪声,加噪。除噪等信号波形

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 6,207
精华内容 2,482
关键字:

matlab对信号滤波的程序

matlab 订阅