分幀和分窗處理：

• 對信號x加窗分幀處理

  1 2 3 4 5 6 7

  wlen=50; % 幀長 inc=18; % 幀移 win=hanning(wlen); % 窗函數 fn=floor(((N-wlen)/inc))+1; % 計算幀數 time=(0:N-1)/Fs; frameTime=(((1:fn)-1)*inc+wlen/2)/Fs; X=enframe(x,win,inc)'; %分窗處理

時域特征：

• 求信號signal的時域特征

  • 時域函數

    1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

    function [ rms,peak,factor] = timechara(signal) % timechara:作用，求輸入信號的自相關序列的相關參數函數 % 輸入：signal：一幀信號的序列 % 輸出：rms：每幀信號的有效值（對時間的均值） % peak：每幀信號的峰值 % factor(1):峰值因子，反映波形的形狀特征 % factor(2)：脈沖因子 % factor(3)：裕度因子 % factor(4)：波形因子 % factor(5)：K因子 % kurtosis:峭度 wlen = length(signal(:,1)); % 每幀信號的采樣點個數 s_ave = mean(signal); % 每幀信號的均值 rms = sqrt(sum((signal-s_ave).^2)/wlen); % 每幀信號的自相關序列的有效值（對時間的均值） peak = (max(signal)-min(signal))/2; % 每幀信號的峰值 factor(1) = peak/rms; % crestfactor：每幀信號的峰值因子，反映波形的形狀特征; factor(2) = peak/(sum(abs(signal))/wlen); % impulsefactor:每幀信號的脈沖因子 factor(3) = peak/((sum(sqrt(abs(signal)))/wlen).^2); % marginfactor:每幀信號的裕度因子 factor(4) = rms/(sum(abs(signal))/wlen); % shapefactor:每幀信號的波形因子 factor(5) = peak*rms; % Kfactor:每幀信號的K因子 factor(6) = kurtosis(signal); % kurtosis:峭度 end

  • 例子

    1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

    x=0:0.1:2*pi; y=sin(x); %信號 ma = max(y); %最大值 mi = min(y); %最小值 me = mean(y); %平均值 pk = ma-mi; %峰-峰值 av = mean(abs(y)) %絕對值的平均值(整流平均值) va = var(y); %方差 st = std(y); %標準差 ku = kurtosis(y); %峭度 rm = rms(y); %均方根 S = rm/av %波形因子 C = pk/rm; %峰值因子 Kr = sum(y.^4)/sqrt(sum(y.^2)) %峭度因子 I = pk/av %脈沖因子 xr = mean(sqrt(abs(y)))^2; L = pk/xr; %裕度因子

• 求信號的短時平均能量的函數。

  1
  2
  3
  4
  5
  6
  7
  8
  function [ E ] = energy( frame )
  % energy:對于一幀信號，求出短時能量
  % 輸入參數：frame:一幀信號
  % 輸出參數：E:信號的短時能量
  u=frame;
  u2=u.*u;
  E=sum(u2); %幅度的平方和
  end
  

• 短時平均過零率

  1
  2
  3
  4
  5
  6
  7
  8
  9
  10
  11
  12
  13
  function [ zcr ] = zcr( frame )
  % zcr:求一幀信號的短時平均過零數
  % 輸入參數：
  zcr = 0;
  wlen = length(frame);
  frame = frame-mean(frame); % 消除直流分量,但是此處處理后與在整個信號出消除直流分量求得的zrc有差距
  for i=1: (wlen-1) % 在一幀內尋找過零點
  %if條件處可加門限
  if frame(i)* frame(i+1)< 0 % 判斷是否為過零點
  zcr = zcr+1; % 是過零點，記錄1次
  end
  end
  end
  

• 求輸入信號的自相關序列的相關參數函數

  1
  2
  3
  4
  5
  6
  7
  8
  9
  10
  11
  12
  13
  14
  15
  function [ rms,peak,crestfactor] = fcorr(signal)
  % fcorr:作用，求輸入信號的自相關序列的相關參數函數
  % 輸入：signal一個信號的序列
  % 輸出：rms：每個信號的自相關序列的有效值（對時間的均值）
  % peak：自相關函數序列的峰值
  % crestfactor：自相關函數序列的峰值因子，反映波形的形狀特征
  wlen = length(signal(:,1));
  ss = xcorr(signal);
  ss = ss(1:wlen,:); % 因為相關函數是對稱的，取前一半
  ss_ave = mean(ss); % 自相關函數序列的均值
  % rms = sqrt(sum((signal-ss_ave).^2)/wlen); % 這是rms的錯誤求法，但是測試的訓練效果較好
  rms = sqrt(sum((ss-ss_ave).^2)/wlen); % 每個信號的自相關序列的有效值（對時間的均值）
  peak = (max(ss)-min(ss))/2; % 峰值
  crestfactor = peak/rms; % 峰值因子，反映波形的形狀特征;
  end
  

時頻域特征

• 求信號的小波的特性

  1 2 3 4 5 6 7 8 9 10 11 12 13

  function [ energy,sqr ] = wave( frame ) % 輸入參數：frame:一幀信號 % 輸出參數： %---------energy:小波系數的能量 %---------sqr:小波系數的均方差 T = wpdec(frame,2,'db2'); %db8小波做兩層分解 %重構第二層的所有節點 for i = 1:4 y = wprcoef(T,[2,i-1]); energy(:,i) = sum(y.*y); %重構信號的能量 sqr(:,i) = var(y); %重構信號的均方差 end end

• 小波變換的改進函數

  1 2 3 4 5 6 7 8 9 10 11 12 13 14

  function [ energy,sqr] = wave( frame,N ) % wave:求信號的小波的特性 % 輸入參數：frame:一幀信號 % 輸出參數：energy:小波系數的能量；sqr:小波系數的均方差;coef:波峰系數 T = wpdec(frame,N,'db8'); %db8小波做N層分解 %重構第N層的指定節點 for i = 1:(2^N) y = wprcoef(T,[N,i-1]); % y = wprcoef(T,[N,mum(i)]); energy(:,i) = sum(y.*y); %重構信號的能量 sqr(:,i) = var(y); %重構信號的均方差 % coef(:,i) = max(y)/mean(y); %重構信號的波峰系數，效果不好 end end

• HHT變換的相關特征求取函數

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62

  % 對一個信號做HHT變換，求得相關的特征 function [peaks,bjpvar,aveA,mseA,rmsA,msize,ratio1,msee] = infreqfeature(x,fs) % 輸入：x 是一個信號，是一個列向量 % 輸出： % peaks 邊際譜的前4個峰值 % bjpvar 邊際譜的方差 % aveA 各分量瞬時幅度的均值 % rmsA 各分量瞬時幅度的有效值 % mseeA 各分量瞬時幅度的方差和 % msize(1)、msize(2)、msize(3)、msize(4):瞬時頻率第三個分量的極大值、極小值點數和第四個分量的極大值、極小值點的個數 % ratio 各分量瞬時頻率的方差貢獻率 % % % 調用的函數： % hht工具箱函數： emd()、hhspectrum()、toimage() % 自定義函數： findmainfreq() N = length(x); imp=emd(x); % 對信號進行EMD分解 [m,n]=size(imp); % 求取EMD分解成幾個分量 % 對IMF分量求取瞬時頻率與振幅：A：是每個IMF的振幅向量,f:每個IMF對應的歸一化瞬時頻率，t：時間序列號 [A,f,t] = hhspectrum(imp); % 求瞬時頻率 infreq = fs * f; %求信號的邊界譜bjp，E：對應的振幅值 [E,tt1]=toimage(A,f,t,length(t)); enery = E; E=flipud(E); for k=1:size(E,1) bjp(k)=sum(E(k,:))*1/fs; end % 求邊界譜的峰值和對應的采樣點的位置 % peaks表示前5個峰值，如果需要求這些峰值對應的位置，直接在peaks參數后面加一個參數即可，詳見該函數的定義 [peaks] = findmainfreq(bjp); bjpvar = var(bjp); %邊際譜的方差 % 瞬時頻率第三個和第四個分量的極值點的個數 [indmin3, indmax3, indzer] = extr(infreq(3,:),1/fs); msize(1) = length(indmin3); % 第三個分量極大值點數 msize(2) = length(indmax3); % 第三個分量極小值點數 [indmin4, indmax4, indzer] = extr(infreq(4,:),1/fs); msize(3) = length(indmin4); % 第四個分量極大值點數 msize(4) = length(indmax4); % 第四個分量極大值點數 % 非有效特征，瞬時頻率的均值、方差和有效值的特征對于好壞瓶子的區分性很差，所以不需要使用 % 所有分量瞬時頻率的均值、方差和有效值 for k = 1:m-1 % ave(k) = mean(infreq(k,:)); % 瞬時頻率的均值 mse(k) = var(infreq(k,:)); % 各分量的瞬時頻率的方差 % rms(k) = sqrt(sum((infreq(k,:) - ave(k)).^2)/N);% 各分量的瞬時頻率的有效值 end % 所有分量瞬時幅度的均值、方差和有效值 for k = 1:m-1 aveA(k) = mean(A(k,:)); % 瞬時幅度的均值 mseA(k) = var(A(k,:)); % 各分量的瞬時幅度的方差 rmsA(k) = sqrt(sum((A(k,:) - aveA(k)).^2)/N); % 各分量的瞬時幅度的有效值 end % 所有分量的方差和 msee = sum(mse); % 各分量的方差貢獻率,此處只求前6個分量的方差貢獻率 for j = 1:6 ratio(j) = (100*mse(j))/msee; %獲取瞬時頻率的方差貢獻率 end ratio1 = (100*mse(1))/msee; % 第一個分量的瞬時頻率方差貢獻率 end

頻域特征

• 求出信號的頻域特征參數

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

  function [ area,peaks,loc ] = spectrum( frame ) % spectrum:對一幀信號取頻譜特征 % 輸入參數：frame:一幀信號，長度為256 % 輸出參數：area:頻譜面積；peaks：頻譜曲線峰值系數(最大5個)：col:最大5個峰值系數對應的位置 N = length(frame); [f amp] = FFT(frame,N); temp = amp(1:N/2)'; area = polyarea(1:N/2,temp); %求頻譜面積 [K,V] = findpeaks(temp); %找出頻譜峰值 peaks = sort(K,'descend'); %降序排列峰值 peaks = peaks(1:5); %取出前5大的峰值 for i = 1:5 mark = find(K == peaks(i)); loc(i) = V(mark); end end

• 直接求頻域特征

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

  function [ FC,MSF,VF,RMSF,RVF ] = frequency( frame ) % frequency:求出一幀信號的頻域特征參數 % 輸入參數：frame:一幀信號 % 輸出參數： FC:重心頻率 % MSF:均方頻率 % VF:頻率方差 % RMSF:均方根頻率 % RVF:頻率標準差 delta = 1/48000; % 采樣間隔 N = length(frame(:,1)); % 每幀信號的點數 for i=2:N s(i,:) = (frame(i,:)-frame(i-1,:))/delta; end FC = sum(s(2:N,:).*frame(2:N,:))/(2*pi*sum(frame(:,:).^2)); MSF = sum(s(2:N,:).^2)/(4*(pi^2)*sum(frame(:,:).^2)); VF = MSF-(FC^2); RMSF = sqrt(MSF); RVF = sqrt(VF); end

參考：音頻特征提取——常用音頻特征

總結

以上是生活随笔為你收集整理的MATLAB求音频信号特征的自定义函数.md的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。