一、簡介

1 心音：
心臟收縮舒張時產生的聲音，可用耳或聽診器在胸壁聽到，亦可用電子儀器記錄下來（心音圖）。可分為第一心音（S1）第二心音（S2）。（正常情況下均可聽到）。第三心音(S3通常僅在兒童及青少年可聽到)，第四心音（S4正常情況很少聽到）。從心臟產生的心音經過組織的介導傳到胸壁表面，其中以骨傳導最好。
心音是心臟及心血管系統機械運動狀況的反映，其中包含著心臟各個部分本身及相互之間作用的生理和病理資訊。心音訊號的識別與分類對心血管系統疾病的診斷具有重要的意義，其準確性、可靠性的好壞決定著診斷與治療心臟病患者的效果。早期的心音識別與分類是醫生根據聽診結果來完成的，顯然這一過程具有一定的主觀性且可靠性不高。隨著訊號處理與分析技術的不斷髮展，對心音的研究也逐步由定性分析進入了定量分析的階段。國內外許多生物醫學工程研究人員將傳統的模式識別方法，以及神經網路方法用於心音的識別與分類，期望實現心音的自動解釋和自動診斷，以便向臨床醫生提供實用的輔助診斷資訊。此外，心音的識別與分類還有助於對心音產生機制的理解。
心音訊號研究主要是採用微電子技術，檢測技術，現代數字訊號處理技術和生物醫學工程技術，研究和揭示心音與心臟病之間的關係。
無論是多大單位的記錄，人的體重，心率，血壓等生理引數，都是時變的，稱為心率或者血壓的變異性。
心音是在心動週期中，由於心肌收縮和舒張，瓣膜啟閉，血流衝擊心室壁和大 動脈等因素引起的機械振動，通過周圍組織傳到胸壁，將耳緊貼胸壁或將聽診器放在胸壁 定部位，聽到的聲音 通常很容易聽到第一和第二心音，有時在某些情況下聽到第三或第四心音。
第一心音：S1,發生在心臟收縮期開始，音調低沉，持續時間較長(約0.15秒)。產生的原因包括心室肌的收縮，房室瓣突然關閉以及隨後射血入主動脈等引起的振動。第一心音的最佳聽診部位在鎖骨中線第五肋間隙或在胸骨右緣。相對於心電圖上QRS波開始後0.02_{0.04s，佔時0.08}

心音雜音對正常心音形成了一定的干擾，但心音雜音的出現對心音訊號分析具有實際應用價值和臨床意義。
根據雜音出現的時間與S1,S2心音的關係，可分為早，中，晚期雜音。雜音的強度一般可視其振幅與S1比較分類。大於S1,強。小於S1,大於1/3S1，中。低於S1的1/3則為低，僅有輕微振動則為極低振幅雜音。
由於心音訊號屬於強噪聲背景下的人體微弱生物訊號，由於心音訊號是由複雜的生命體發出的不穩定的自然訊號。心音的改變和雜音的出現往往是心臟產生器質性病變的早期症狀，心臟內部的物理結構發生變化將直接影響和改變心音訊號。

目前廣泛應用 的心電圖檢查是心臟變時性和變傳導性的最佳檢測方法，但不能用來檢測心臟的變力性先天心臟瓣膜受損。心臟傳導組織病變引起的心臟機械活動障礙不會首先反應在心電圖上，卻能首先反應在心音訊號上。當冠心病的阻塞達到70%以上時才能引起心電圖訊號的改變，實際上，達到25%就可以改變心音訊號。

2 異常心音
　　包括S2、S2的異常及收縮期，舒張期的附加音（或額外音）。
　　① 第一心音異常。指S1增強，減弱或分裂。估計S1的響度，最好是用聽診，心音圖判斷能力有限。S1增強、減弱或強弱不等的臨床情況見表1。S1分裂指M1與T1相距>0.04sec，可見於正常兒童、青年及體瘦者，無重要意義。S1異常分裂時M1與T1相距可大於0.06秒，見於電激動延遲(如右束支傳導阻滯等)、機械活動延遲（如房間隔缺損、嚴重二尖瓣狹窄等。聽診S1分裂在二尖瓣及三尖瓣區最清楚坐位及呼氣時更清楚。
　　② 第二心音異常。包括S2增強，減弱或分裂。
　　S2增強又分P2增強及A2增強，P2增強見於肺血流量增多（如間隔缺損），肺血管阻力增加；肺靜脈壓力增高(如二尖瓣狹窄)。P2亢進一般在肺動脈瓣區聽到。A2增強見於體迴圈阻力增高或血流量增多，向肺動脈瓣及心尖區傳導，見於 高血壓等。
　　S2減弱又分P2減弱及A2減弱。P2減弱見於肺動脈壓低、肺血流量減少或 肺動脈瓣狹窄等。A2減弱見於體迴圈阻力低、血流減少、血壓低、 主動脈瓣狹窄或嚴重關閉不全。
　　S2分裂可為生理性。右室射血結束稍晚於左室，吸氣時P2延遲出現，此時可聽到S2分裂。但呼氣時A2與P2接近或重疊，分裂消失。這見於青少年，在肺動脈瓣區聽診明顯，坐位時可消失。
　　S2異常分裂包括寬分裂、固定分裂、逆分裂、分裂減窄。寬分裂是呼氣時S2分裂不消失，見於右室射血時間延長或左室射血時間縮短等情況。S2固定分裂指呼吸時A2 -P2間期無明顯改變或變動<0.02sec，見於分流量較大的房間隔缺損等。S2逆分裂指A2在P2之後，吸氣時A2-P2分裂不顯，呼氣時P2提早出現，分裂增寬。見於主動脈瓣關閉延遲。S2分裂減窄常由於嚴重肺動脈高壓P2較早出現所致。

3 心音訊號分析方法：
傳統的譜分析方法通過快速傅立葉變換將時，頻域關聯起來。但FFT時頻域分離，並以訊號的頻率特性時不變，或統計特性穩定為前提。
傳統的穩態分析方法反映的是訊號的靜態頻譜特徵，對於包括人體心音訊號在內的生物學生理訊號，由於環境的影響而表現為非平穩時變特性。因此採用經典譜分析方法難以準確反映心音訊號的動態變化。
在傳統心音分析的基礎上，提出了很多方法：
1.短時傅立葉變換（STFT）
2.小波變換
3.其他時頻分析方法

二、原始碼

clc;
clear;
close all;

%% 提取訊號
M = importdata('3.txt');
fsample=1000;%取樣率為1KHz

[mx,my]=size(M);
Signal=M(:,2);%M的第一列為時間，第二列為訊號

length=floor(mx/2);%取原始訊號的一半。
S=Signal(1:length);

%% 高通濾波，去除基線漂移的影響
disp('-------------------------------------------');
disp('1：工具箱巴特沃斯高通濾波器');
disp('2：IIR高通濾波');
disp('3：FIR高通濾波');
disp('4：中值濾波');
disp('5：稀疏小波濾波');
disp('6：中值+小波濾波');
disp('-------------------------------------------');
choose=input('選擇濾波方式choose=');
function [s, err_mse, iter_time]=sparseOMP(x,A,m,varargin)
[n1 n2]=size(x);
if n2 == 1
    n=n1;
elseif n1 == 1
    x=x';
    n=n2;
else
   error('x must be a vector.');
end
    
sigsize     = x'*x/n;
initial_given=0;  
err_mse     = [];
iter_time   = [];
STOPCRIT    = 'M';
STOPTOL     = ceil(n/4);
MAXITER     = n;
verbose     = false;
s_initial   = zeros(m,1);
GradSteps   = 'auto';
alpha       = 1;
weakness    = 1;

if verbose
   display('Initialising...') 
end

switch nargout 
    case 3
        comp_err=true;
        comp_time=true;
    case 2 
        comp_err=true;
        comp_time=false;
    case 1
        comp_err=false;
        comp_time=false;
    case 0
        error('Please assign output variable.')        
    otherwise
        error('Too many output arguments specified')
end


% Put option into nice format
Options={};
OS=nargin-3;
c=1;
for i=1:OS
    if isa(varargin{i},'cell')
        CellSize=length(varargin{i});
        ThisCell=varargin{i};
        for j=1:CellSize
            Options{c}=ThisCell{j};
            c=c+1;
        end
    else
        Options{c}=varargin{i};
        c=c+1;
    end
end
OS=length(Options);
if rem(OS,2)
   error('Something is wrong with argument name and argument value pairs.') 
end
for i=1:2:OS
   switch Options{i}
        case {'stopCrit'}
            if (strmatch(Options{i+1},{'M'; 'corr'; 'mse'; 'mse_change'},'exact'));
                STOPCRIT    = Options{i+1};  
            else error('stopCrit must be char string [M, corr, mse, mse_change]. Exiting.'); end 
        case {'stopTol'}
            if isa(Options{i+1},'numeric') ; STOPTOL     = Options{i+1};   
            else error('stopTol must be number. Exiting.'); end
        case {'P_trans'} 
            if isa(Options{i+1},'function_handle'); Pt = Options{i+1};   
            else error('P_trans must be function _handle. Exiting.'); end
        case {'maxIter'}
            if isa(Options{i+1},'numeric'); MAXITER     = Options{i+1};             
            else error('maxIter must be a number. Exiting.'); end
        case {'wf'}
            if isa(Options{i+1},'numeric'); alpha       = Options{i+1}; 
                if alpha <1 weakness =0; else alpha =1; weakness = 1; end          
            else error('wf must be a number. Exiting.'); end
        case {'verbose'}
            if isa(Options{i+1},'logical'); verbose     = Options{i+1};   
            else error('verbose must be a logical. Exiting.'); end 
        case {'start_val'}
            if isa(Options{i+1},'numeric') && length(Options{i+1}) == m ;
                s_initial     = Options{i+1};   
                initial_given=1;
            else error('start_val must be a vector of length m. Exiting.'); end
        case {'GradSteps'}
            if isa(Options{i+1},'numeric') || strcmp(Options{i+1},'auto') ;
                GradSteps     = Options{i+1};   
            else error('start_val must be a vector of length m. Exiting.'); end
        otherwise
            error('Unrecognised option. Exiting.') 
   end
end




if strcmp(STOPCRIT,'M') 
    maxM=STOPTOL;
else
    maxM=MAXITER;
end

if nargout >=2
    err_mse = zeros(maxM,1);
end
if nargout ==3
    iter_time = zeros(maxM,1);
end


if          isa(A,'float')      P =@(z) A*z;  Pt =@(z) A'*z;
elseif      isobject(A)         P =@(z) A*z;  Pt =@(z) A'*z;
elseif      isa(A,'function_handle') 
    try
        if          isa(Pt,'function_handle'); P=A;
        else        error('If P is a function handle, Pt also needs to be a function handle. Exiting.'); end
    catch error('If P is a function handle, Pt needs to be specified. Exiting.'); end
else        error('P is of unsupported type. Use matrix, function_handle or object. Exiting.'); end

if initial_given ==1;
    IN          = find(s_initial);
    Residual    = x-P(s_initial);
    s           = s_initial;
    oldERR      = Residual'*Residual/n;
else
    IN          = [];
    Residual    = x;
    s           = s_initial;
    sigsize     = x'*x/n;
    oldERR      = sigsize;
end

        mask=zeros(m,1);
        mask(ceil(rand*m))=1;
        nP=norm(P(mask));
        if abs(1-nP)>1e-3;
            display('Dictionary appears not to have unit norm columns.')
        end

if verbose
   display('Main iterations...') 
end
tic
t=0;
normA=(sum(A.^2,1)).^0.5;
NI = 1:size(A,2);
p=zeros(m,1);
DR=Pt(Residual);
[v I]=max(abs(DR(NI))./normA(NI)');
INI = NI(I);
IN=[IN INI];
NI(I) = [];
if weakness ~= 1
    [vals inds] = sort(abs(DR),'descend');
    I=inds( find( vals >= alpha * v ) );
end
    
IN = union(IN,I);
if strcmp(STOPCRIT,'M') & length(IN) >= STOPTOL
    IN=IN(1:STOPTOL);
end
MASK=zeros(size(DR));
pDDp=1;
done = 0;
iter=1;

while ~done

    
    % Select new element
    if isa(GradSteps,'char')
        if strcmp(GradSteps,'auto')
             
%             finished=0;    
%             while ~finished
            % Update direction    
                 if iter==1
                     p(IN)=DR(IN);
                     Dp=P(p);
                 else
                     MASK(IN)=1;
                     PDR=P(DR.*MASK);
                     b=-Dp'*PDR/pDDp;
                     p(IN)=DR(IN) +b*p(IN);
                     Dp=PDR +b* Dp;
                 end
             % Step size
%                  Dp=P(p); % =P(DR(IN)) +b P(p(IN));
                 pDDp=Dp'*Dp;
                 a=Residual'*Dp/(pDDp);
             % Update coefficients   
                 s=s+a*p;
             % New Residual and inner products
                 Residual=Residual-a*Dp;
                 DR=Pt(Residual);
                 % select new element
                     [v I]=max(abs(DR(NI))./normA(NI)');
                     INI = NI(I);
                     if weakness ~= 1
                         [vals inds] = sort(abs(DR),'descend');
                         I=inds( find( vals >= alpha * v ) );
                     end
                     IN = union(IN,INI);
                     NI(I) = [];
                     if strcmp(STOPCRIT,'M') & length(IN) >= STOPTOL
                        IN=IN(1:STOPTOL);
                     end


        else
            
            
            error('Undefined option for GradSteps, use ''auto'' or an integer.')
        end
    elseif isa(GradSteps,'numeric') 

                
        % Do GradSteps gradient steps
        count=1;
        while count<=GradSteps
            % Update direction    
                 if iter==1
                     p(IN)=DR(IN);
                     Dp=P(p);
                 else
                     MASK(IN)=1;
                     PDR=P(DR.*MASK);
                     b=-Dp'*PDR/pDDp;
                     p(IN)=DR(IN) +b*p(IN);
                     Dp=PDR +b* Dp;
                 end
             % Step size
%                  Dp=P(p);   
                 pDDp=Dp'*Dp;
                 a=Residual'*Dp/(pDDp);
             % Update coefficients   
                 s=s+a*p;
             % New Residual and inner products
                 Residual=Residual-a*Dp;
                 DR=Pt(Residual);
                  count=count+1;
        end
             % select new element
                 [v I]=max(abs(DR(NI))./normA(NI)');
                 INI = NI(I);
                 if weakness ~= 1
                     [vals inds] = sort(abs(DR),'descend');
                     I=inds( find( vals >= alpha * v ) );
                 end
                 IN = union(IN,INI);
                  NI(I) = [];
                 if strcmp(STOPCRIT,'M') & length(IN) >= STOPTOL
                    IN=IN(1:STOPTOL);
                 end
                
     else
          error('Undefined option for GradSteps, use ''auto'' or an integer.')
     end


     ERR=Residual'*Residual/n;
     if comp_err
         err_mse(iter)=ERR;
     end
     
     if comp_time
         iter_time(iter)=toc;
     end


     if strcmp(STOPCRIT,'M')
         if length(IN) >= STOPTOL
             done =1;
         elseif verbose && toc-t>10
            display(sprintf('Iteration %i. --- %i selected elements',iter ,length(IN))) 
            t=toc;
         end
    elseif strcmp(STOPCRIT,'mse')
         if comp_err
            if err_mse(iter)<STOPTOL;
                done = 1; 
            elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse',iter ,err_mse(iter))) 
                t=toc;
            end
         else
             if ERR<STOPTOL;
                done = 1; 
             elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse',iter ,ERR)) 
                t=toc;
             end
         end
     elseif strcmp(STOPCRIT,'mse_change') && iter >=2
         if comp_err && iter >=2
              if ((err_mse(iter-1)-err_mse(iter))/sigsize <STOPTOL);
                done = 1; 
             elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse change',iter ,(err_mse(iter-1)-err_mse(iter))/sigsize )) 
                t=toc;
             end
         else
             if ((oldERR - ERR)/sigsize < STOPTOL);
                done = 1; 
             elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse change',iter ,(oldERR - ERR)/sigsize)) 
                t=toc;
             end
         end
     elseif strcmp(STOPCRIT,'corr') 
          if max(abs(DR)) < STOPTOL;
             done = 1; 
          elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i corr',iter ,max(abs(DR)))) 
                t=toc;
          end
     end
     
    % Also stop if residual gets too small or maxIter reached
     if comp_err
         if err_mse(iter)<1e-16
             display('Stopping. Exact signal representation found!')
             done=1;
         end
     else


         if iter>1
             if ERR<1e-16
                 display('Stopping. Exact signal representation found!')
                 done=1;
             end
         end
     end

     if iter >= MAXITER
         display('Stopping. Maximum number of iterations reached!')
         done = 1; 
     end
     
   
     if ~done
        iter=iter+1;
        oldERR=ERR;
     end
end



if nargout >=2
    err_mse = err_mse(1:iter);
end
if nargout ==3
    iter_time = iter_time(1:iter);
end
if verbose
   display('Done') 
end

三、執行結果

四、備註

版本：2014a
完整程式碼或代寫加1564658423