Undefined function 'mcra2_estimation' for input arguments of type 'struct'.

function parameters = initialise_parameters(ns_ps,Srate,method) len_val = length(ns_ps); switch lower(method) case 'martin' L_val=len_val; R_val=len_val/2; D_val=150; V_val=15; Um_val=10; Av_val=2.12; alpha_max_val=0.96; alpha_min_val=0.3; beta_max_val=0.8; x_val=[1 2 5 8 10 15 20 30 40 60 80 120 140 160]; Y_M_val=[0 .26 .48 .58 .61 .668 .705 .762 .8 .841 .865 .89 .9 .91]; Y_H_val=[0 .15 .48 .78 .98 1.55 2.0 2.3 2.52 2.9 3.25 4.0 4.1 4.1]; xi_val=D_val; M_D_val=interp1(x_val,Y_M_val,xi_val); H_D_val=interp1(x_val,Y_H_val,xi_val); xi_val=V_val; M_V_val=interp1(x_val,Y_M_val,xi_val); H_V_val=interp1(x_val,Y_H_val,xi_val); minact_val(1:L_val,1:Um_val)=max(ns_ps); parameters = struct('n',2,'len',len_val,'alpha_corr',0.96,'alpha',0.96*ones(len_val,1),'P',ns_ps,'noise_ps',ns_ps,'Pbar',ns_ps,... 'Psqbar',ns_ps,'actmin',ns_ps,'actmin_sub',ns_ps,'Pmin_u',ns_ps,'subwc',2,'u',1,'minact',minact_val,'lmin_flag',zeros(len_val,1),... 'L',L_val,'R',R_val,'D',D_val,'V',V_val,'Um',Um_val,'Av',Av_val,'alpha_max',alpha_max_val,'alpha_min',alpha_min_val,... 'beta_max',beta_max_val,'Y_M',Y_M_val,'Y_H',Y_H_val,'M_D',M_D_val,'H_D',H_D_val,'M_V',M_V_val,'H_V',H_V_val); case 'mcra' parameters = struct('n',2,'len',len_val,'P',ns_ps,'Pmin',ns_ps,'Ptmp',ns_ps,'pk',zeros(len_val,1),'noise_ps',ns_ps,... 'ad',0.95,'as',0.8,'L',round(1000*2/20),'delta',5,'ap',0.2); case 'imcra' alpha_d_val=0.85; alpha_s_val=0.9; U_val=8;V_val=15; Bmin_val=1.66;gamma0_val=4.6;gamma1_val=3; psi0_val=1.67;alpha_val=0.92;beta_val=1.47; j_val=0; b_val=hanning(3); B_val=sum(b_val); b_val=b_val/B_val; Sf_val=zeros(len_val,1);Sf_tild_val=zeros(len_val,1); Sf_val(1) = ns_ps(1); for f=2:len_val-1 Sf_val(f)=sum(b_val.*[ns_ps(f-1);ns_ps(f);ns_ps(f+1)]); end Sf_val(len_val)=ns_ps(len_val); Sf_tild_val = zeros(len_val,1); parameters = struct('n',2,'len',len_val,'noise_ps',ns_ps,'noise_tild',ns_ps,'gamma',ones(len_val,1),'Sf',Sf_val,... 'Smin',Sf_val,'S',Sf_val,'S_tild',Sf_val,'GH1',ones(len_val,1),'Smin_tild',Sf_val,'Smin_sw',Sf_val,'Smin_sw_tild',Sf_val,... 'stored_min',max(ns_ps)*ones(len_val,U_val),'stored_min_tild',max(ns_ps)*ones(len_val,U_val),'u1',1,'u2',1,'j',2,... 'alpha_d',0.85,'alpha_s',0.9,'U',8,'V',15,'Bmin',1.66,'gamma0',4.6,'gamma1',3,'psi0',1.67,'alpha',0.92,'beta',1.47,... 'b',b_val,'Sf_tild',Sf_tild_val); case 'doblinger' parameters = struct('n',2,'len',len_val,'alpha',0.7,'beta',0.96,'gamma',0.998,'noise_ps',ns_ps,'pxk_old',ns_ps,... 'pxk',ns_ps,'pnk_old',ns_ps,'pnk',ns_ps); case 'hirsch' parameters = struct('n',2,'len',len_val,'as',0.85,'beta',1.5,'omin',1.5,'noise_ps',ns_ps,'P',ns_ps); case 'mcra2' freq_res=Srate/len_val; k_1khz=floor(1000/freq_res); k_3khz=floor(3000/freq_res); delta_val=[2*ones(k_1khz,1);2*ones(k_3khz-k_1khz,1);5*ones(len_val/2-k_3khz,1);... 5*ones(len_val/2-k_3khz,1);2*ones(k_3khz-k_1khz,1);2*ones(k_1khz,1)];

        parameters = struct('n',2,'len',len_val,'ad',0.95,'as',0.8,'ap',0.2,'beta',0.8,'beta1',0.98,'gamma',0.998,'alpha',0.7,...
            'delta',delta_val,'pk',zeros(len_val,1),'noise_ps',ns_ps,'pxk_old',ns_ps,'pxk',ns_ps,'pnk_old',ns_ps,'pnk',ns_ps);

      case 'conn_freq'
        D = 7; 
        b = triang(2*D+1)/sum(triang(2*D+1));
        b = b';
        beta_min = 0.7; % for R's recursion
        U = 5;
        V = 8;
        gamma1 = 6; 
        gamma2 = 0.5; 
        K_tild = 2*sum(b.^2)^2/sum(b.^4);
        alpha_max_val=0.96; 
        alpha_min_val=0.3;
        stored_min = max(ns_ps)*ones(len_val,U);

        alpha_c = 0.7;
        noise_ps = ns_ps;
        Rmin_old = 1;
        Pmin_sw = ns_ps;
        Pmin = ns_ps;
        P = ns_ps;
        Decision = zeros(size(P));
        u1 = 1;
        j = 0;
        parameters = struct('len',len_val,'D',D,'b',b,'U',U,'V',V,'gamma1',gamma1,'gamma2',gamma2,'K_tild',K_tild,'alpha_c',alpha_c,...
            'noise_ps',noise_ps,'Rmin_old',Rmin_old,'Pmin_sw',Pmin_sw,'Pmin',Pmin,'SmthdP',P,'u1',u1,'j',j,'alpha',0,'alpha_max',alpha_max_val,...
            'stored_min',stored_min,'beta_min',beta_min,'Decision',Decision);

    otherwise
            error('Method not implemented. Check spelling.');
end
function [parameters]=martin_estimation(ns_ps,parameters);

YFRAME = ns_ps; alpha_corr = parameters.alpha_corr; alpha = parameters.alpha; P = parameters.P; Pbar = parameters.Pbar; Psqbar = parameters.Psqbar; actmin = parameters.actmin; actmin_sub = parameters.actmin_sub; minact = parameters.minact; Pmin_u = parameters.Pmin_u; subwc = parameters.subwc; u = parameters.u; lmin_flag = parameters.lmin_flag; n = parameters.n; L = parameters.L; R = parameters.R; Um = parameters.Um; V = parameters.V; D = parameters.D; Av = parameters.Av; alpha_max = parameters.alpha_max; alpha_min = parameters.alpha_min; beta_max = parameters.beta_max; M_D = parameters.M_D; M_V = parameters.M_V; H_D = parameters.H_D; H_V = parameters.H_V; noise_est = parameters.noise_ps;

%calculating the optimal smoothing correction factor alpha_corr_t=1/(1+ ((sum(P)/sum(YFRAME)-1)^2)); alpha_corr=0.7*alpha_corr+0.3*max(alpha_corr_t,0.7);

%calculating the optimal smoothing factor

alpha=(alpha_max*alpha_corr)./((P./(noise_est+eps) -1).^2 +1); %

%calculation of smoothed periodogram P=alpha.*P+((1-alpha).*YFRAME); %calculation of variance of P(i,k) and Qeq(i,k) bet=alpha.^2; bet=min(bet,beta_max); Pbar=bet.*Pbar+(1-bet).*P; Psqbar=bet.*Psqbar+(1-bet).*(P.^2); varcap_P=abs(Psqbar-(Pbar.^2)); Qeqinv=varcap_P./(2*(noise_est.^2)); Qeqinv=min(Qeqinv,0.5); Qeq=1./(Qeqinv+eps); %calculation of Bmin(i,k) and Bmin_sub(i,k)

Qeq_tild=(Qeq-2*M_D)./(1-M_D); Qeq_tild_sub=(Qeq-2*M_V)./(1-M_V); %Bmin=1+(((D-1)*2./Qeq_tild).*gamma((1+(2./Qeq)).^H_D)); %Bmin_sub=1+(((V-1)*2./Qeq_tild).*gamma((1+(2./Qeq)).^H_V)); Bmin=1+(D-1)*2./Qeq_tild; % using the approximation in Eq. 17 Bmin_sub=1+(V-1)*2./Qeq_tild_sub;

%calculation of Bc(i) Qinv_bar=(1/L)*sum((1/Qeq)); Bc=1+Av*sqrt(Qinv_bar);

%calculation of actmin(i,k) and actmin_sub(i,k) k_mod=zeros([L 1]); k_mod(find(P.*Bmin.*Bc<actmin))=1; actmin_sub(find(k_mod))=P(find(k_mod)).*Bmin_sub(find(k_mod)).*Bc; actmin(find(k_mod))=P(find(k_mod)).*Bmin(find(k_mod)).*Bc;

if subwc==V

    %check whether the minimum is local minimum or not
    lmin_flag(find(k_mod))=0;

    %storing the value of actmin(i,k)
    minact(:,u)=actmin;

    %calculation of Pmin_u the minimum of the last U stored values of actmin
    Pmin_u=min(minact,[],2);

    %calculation of noise slope max
    if Qinv_bar<0.03
        noise_slope_max=8;
    elseif Qinv_bar<0.05
        noise_slope_max=4;
    elseif Qinv_bar<0.06
        noise_slope_max=2;
    else
        noise_slope_max=1.2;
    end

    %update Pmin_u if the minimum falls within the search range
    test=find((lmin_flag & (actmin_sub<(noise_slope_max*Pmin_u)) & (actmin_sub>Pmin_u)));
    Pmin_u(test)=actmin_sub(test);
    %noise_est=min(actmin_sub,Pmin_u);
    for x=1:Um
        minact(test,x)=actmin_sub(test);
    end
    actmin(test) = actmin_sub(test);
    lmin_flag(:)=0;
    subwc=1;
    actmin=P;
    actmin_sub=P;   
    if u == Um
        u=1;
    else
        u=u+1;
    end
else
    if subwc>1
        lmin_flag(find(k_mod))=1;
        noise_est=min(actmin_sub,Pmin_u);
        Pmin_u=noise_est;                        
    end
    %noise_est=min(actmin_sub,Pmin_u);    
    subwc=subwc+1;
end

parameters.alpha_corr = alpha_corr; parameters.alpha = alpha; parameters.P = P; parameters.Pbar = Pbar; parameters.Psqbar = Psqbar; parameters.actmin = actmin; parameters.actmin_sub = actmin_sub; parameters.minact = minact; parameters.Pmin_u = Pmin_u; parameters.subwc = subwc; parameters.u = u; parameters.lmin_flag = lmin_flag; parameters.n = n+1; parameters.L = L; parameters.R = R; parameters.Um = Um; parameters.V = V; parameters.D = D; parameters.Av = Av; parameters.alpha_max = alpha_max; parameters.alpha_min = alpha_min; parameters.beta_max = beta_max; parameters.M_D = M_D; parameters.M_V = M_V; parameters.H_D = H_D; parameters.H_V = H_V; parameters.noise_ps = noise_est; function parameters = noise_estimation(noisy_ps,method,parameters) switch lower(method) case 'martin' parameters = martin_estimation(noisy_ps,parameters); case 'mcra' parameters = mcra_estimation(noisy_ps,parameters); case 'imcra' parameters = imcra_estimation(noisy_ps,parameters); case 'doblinger' parameters = doblinger_estimation(noisy_ps,parameters); case 'hirsch' parameters = hirsch_estimation(noisy_ps,parameters); case 'mcra2' parameters = mcra2_estimation(noisy_ps,parameters); case 'conn_freq' parameters = connfreq_estimation(noisy_ps,parameters); end return;function specsub_ns(filename,method,outfile) % % % Implements the basic spectral subtraction algorithm [1] with different % noise estimation algorithms specified by 'method'. % % % Usage: specsub_ns(noisyFile, method, outputFile) % % infile - noisy speech file in .wav format % outputFile - enhanced output file in .wav format % method - noise estimation algorithm: % 'martin' = Martin''s minimum tracking algorithm [2] % 'mcra' = Minimum controlled recursive average algorithm (Cohen) [3] % 'mcra2' = variant of Minimum controlled recursive % average algorithm (Rangachari & Loizou) [4] % 'imcra' = improved Minimum controlled recursive % average algorithm (Cohen) [5] % 'doblinger' = continuous spectral minimum tracking(Doblinger) [6] % 'hirsch' = weighted spectral average (Hirsch & % Ehrilcher) [7] % 'conn_freq' = connected frequency regions (Sorensen & % Andersen) [8] % % % Example call: specsub_ns('sp04_babble_sn10.wav','mcra2','out_ss_mcra2.wav'); % % References: % [1] Berouti, M., Schwartz, M., and Makhoul, J. (1979). Enhancement of speech % corrupted by acoustic noise. Proc. IEEE Int. Conf. Acoust., Speech, % Signal Processing, 208-211. % [2] Martin, R. (2001). Noise power spectral density estimation based on optimal % smoothing and minimum statistics. IEEE Transactions on Speech and Audio % Processing, 9(5), 504-512. % [3] Cohen, I. (2002). Noise estimation by minima controlled recursive averaging % for robust speech enhancement. IEEE Signal Processing Letters, % 9(1), 12-15. % [4] Rangachari, S. and Loizou, P. (2006). A noise estimation algorithm for % highly nonstationary environments. Speech Communication, 28, % 220-231. % [5] Cohen, I. (2003). Noise spectrum estimation in adverse environments: % Improved minima controlled recursive averaging. IEEE Transactions on Speech % and Audio Processing, 11(5), 466-475. % [6] Doblinger, G. (1995). Computationally efficient speech enhancement by % spectral minima tracking in subbands. Proc. Eurospeech, 2, 1513-1516. % [7] Hirsch, H. and Ehrlicher, C. (1995). Noise estimation techniques for robust % speech recognition. Proc. IEEE Int. Conf. Acoust. , Speech, Signal % Processing, 153-156. % [8] Sorensen, K. and Andersen, S. (2005). Speech enhancement with natural % sounding residual noise based on connected time-frequency speech presence % regions. EURASIP J. Appl. Signal Process., 18, 2954-2964. % % Authors: Sundar Rangachari, Yang Lu and Philipos C. Loizou % % Copyright (c) 2006 by Philipos C. Loizou % $Revision: 0.0 $ $Date: 10/09/2006 $ %-------------------------------------------------------------------------

if nargin<3 fprintf('Usage: specsub_ns(infile.wav,method,outfile.wav) \n'); % fprintf(' where ''method'' is one of the following:\n'); % fprintf(' martin = Martin''s minimum tracking algorithm\n'); % fprintf(' mcra = Minimum controlled recursive average algorithm (Cohen) \n') % fprintf(' mcra2 = variant of Minimum controlled recursive average algorithm (Rangachari & Loizou)\n') % fprintf(' imcra = improved Minimum controlled recursive average algorithm (Cohen) \n'); % fprintf(' doblinger = continuous spectral minimum tracking (Doblinger) \n'); % fprintf(' hirsch = weighted spectral average (Hirsch & Ehrilcher) )\n'); % fprintf(' conn_freq = connected frequency regions (Sorensen & Andersen)\n'); % fprintf('\n For more help, type: help specsub_ns\n\n'); return; end [FileName,PathName] = uigetfile('sp04_babble_sn10.wav'); PathOriginal = fullfile(PathName, FileName); [x,Srate,nbits] = wavread(PathOriginal); % PathOriginal = fullfile('C:\Users','sp04_babble_sn10.wav'); % [x,Srate,nbits]=wavread('C:\Users','sp04_babble_sn10.wav');

% =============== Initialize variables =============== %

len=floor(20*Srate/1000); % Frame size in samples if rem(len,2)==1, len=len+1; end; PERC=50; % window overlap in percent of frame size len1=floor(len*PERC/100); len2=len-len1;

alpha=2.0; % power exponent FLOOR=0.002; win=hamming(len); %tukey(len,PERC); % define window

%--- allocate memory and initialize various variables

k=1; nFFT=2*len; img=sqrt(-1); x_old=zeros(len1,1); Nframes=floor(length(x)/len2)-1; xfinal=zeros(Nframes*len2,1);

%=============================== Start Processing ======================================================= % for n=1:Nframes

   insign=win.*x(k:k+len-1);     %Windowing  
   spec=fft(insign,nFFT);     %compute fourier transform of a frame
   sig=abs(spec); % compute the magnitude
   ns_ps=sig.^2;

   % ----------------- estimate/update noise psd --------------
   if n == 1
         parameters = initialise_parameters(ns_ps,Srate,method);   
    else
        parameters = noise_estimation(ns_ps,method,parameters);
   end

   noise_ps = parameters.noise_ps;
   noise_mu=sqrt(noise_ps);  % magnitude spectrum
   % ---------------------------------------------------------

    %save the phase information for each frame.
    theta=angle(spec);  

   SNRseg=10*log10(norm(sig,2)^2/norm(noise_mu,2)^2);

   if alpha==1.0
      beta=berouti1(SNRseg);
   else
     beta=berouti(SNRseg);
  end

   %&&&&&&&&&
   sub_speech=sig.^alpha - beta*noise_mu.^alpha;
   diffw = sub_speech-FLOOR*noise_mu.^alpha;

   % Floor negative components
   z=find(diffw <0);  
   if~isempty(z)
      sub_speech(z)=FLOOR*noise_mu(z).^alpha;
   end

   sub_speech(nFFT/2+2:nFFT)=flipud(sub_speech(2:nFFT/2));  % to ensure conjugate symmetry for real reconstruction
   %multiply the whole frame fft with the phase information
   x_phase=(sub_speech.^(1/alpha)).*(cos(theta)+img*(sin(theta)));

   % take the IFFT 
   xi=real(ifft(x_phase));         

% --- Overlap and add ---------------
% 
xfinal(k:k+len2-1)=x_old+xi(1:len1);
x_old=xi(1+len1:len);

 k=k+len2;
end
%========================================================================================

wavwrite(xfinal,Srate,16,outfile);

%-------------------------------- E N D -------------------------------------- function a=berouti1(SNR)

if SNR>=-5.0 & SNR<=20 a=3-SNR*2/20; else

if SNR<-5.0
 a=4;
end

if SNR>20
  a=1;
end

end

function a=berouti(SNR)

if SNR>=-5.0 & SNR<=20 a=4-SNR*3/20; else

if SNR<-5.0
 a=5;
end

if SNR>20
  a=1;
end

end