EMD和BEMD演算法實現（matlab代碼）

今天結合自己在網上找的一些實現代碼，稍微修改進行測試，沒有進行更多的實驗，可能在一些問題的處理上還是比較毛糙的。
EMD

function imf = emd(x)% Empiricial Mode Decomposition (Hilbert-Huang Transform)% EMD分解或HHT變換% 傳回值為cell類型，依次為一次IMF、二次IMF、...、最後殘差x   = transpose(x(:));imf = [];while ~ismonotonic(x)    x1 = x;    sd = Inf;    while (sd > 0.1) || ~isimf(x1)        s1 = getspline(x1);         % 極大值點樣條曲線        s2 = -getspline(-x1);       % 極小值點樣條曲線        h = x1-(s1+s2)/2;        sd = sum((x1-h).^2)/sum(x1.^2);        x1 = h;    end    imf{end+1} = x1;    x          = x-x1;endimf{end+1} = x;end% 是否單調function u = ismonotonic(x)u1 = length(findpeaks1(x))*length(findpeaks1(-x));if u1 > 0    u = 0;else    u = 1;endend% 是否IMF分量function u = isimf(x)N  = length(x);u1 = sum(x(1:N-1).*x(2:N) < 0);                     % 過零點的個數u2 = length(findpeaks1(x))+length(findpeaks1(-x));    % 極值點的個數if abs(u1-u2) > 1    u = 0;else    u = 1;endend% 據極大值點構造樣條曲線function s = getspline(x)N = length(x);p = findpeaks1(x);s = spline([0 p N+1],[0 x(p) 0],1:N);endfunction n = findpeaks1(x)% Find peaks. 找極大值點，返回對應極大值點的座標n    = find(diff(diff(x) > 0) < 0); % 相當於找二階導小於0的點u    = find(x(n+1) > x(n));n(u) = n(u)+1;  end

BEMD

function [imf_matrix]=bemd(img)%%輸入一副灰階映像[row,col,dep] = size(img);%  row, col and depth of original imageif dep ~= 1    img = im2double(rgb2gray(img));else    img = im2double(img);  end%%%%%主函數% 分解IMF個數設定為3(加上殘餘量為4個分解量)（可根據實際情況修改）m=4;k=1;input_img=img;while(k<m)    [imf_de res_de]=decompsition(input_img);  %% 通過分解得到IMF分量和餘項    imf_matrix(:,:,k)=imf_de;%%儲存IMF分量    input_img=res_de; %%將餘項作為新訊號，再次分解    k=k+1;endimf_matrix(:,:,k)=res_de;%%儲存殘餘量end function  [imf_de res_de]=decompsition(input_img)[width height]=size(input_img);x=1:width;y=1:height;input_img_temple=input_img;while(1)    [zmax imax zmin imin]=extrema2(input_img_temple);  %%%%映像表面極值點    [xmax ymax]=ind2sub(size(input_img_temple),imax);    [xmin ymin]=ind2sub(size(input_img_temple),imin);    [zmaxgrid,~,~]=gridfit(ymax,xmax,zmax,y,x);  %%%%曲面擬合，尋找包絡面的的極值點    [zmingrid,~,~]=gridfit(ymin,xmin,zmin,y,x);    zavggrid=(zmaxgrid+zmingrid)/2;   %%%%包絡均值    %%%%%%IMF分量判斷%%%%%    imf_de=input_img_temple-zavggrid;    SD=sum(sum(imf_de-input_img_temple).^2)/sum(sum(imf_de).^2);    if SD<0.2        break    else         input_img_temple=imf_de;    endendres_de=input_img-imf_de;end

備忘：BEMD利用extrema2尋求曲面極值和gridfit曲面擬合函數實現包絡面的擷取。