acweight_time_filter.m

我认为很不错的语音处理的ｍａｔｌａｂ源代码

function [yAC, errors]=ACweight_time_filter(type, y, Fs, settling_time)
% % ACweight_time_filter: Applies an A or C weighting time filter to a sound recording
% % 
% % Syntax;   
% % 
% % [yAC, errors]=ACweight_time_filter(type, y, Fs, settling_time);
% % 
% % ***********************************************************
% % 
% % Description
% % 
% % This program applies an A or C-weighting filter to 
% % a sound pressure time record, then it outputs the A or C-weighted tme
% % record respectively.
% % 
% % ACweight_time_filter uses the signal processing toolkit.
% % 
% % ***********************************************************
% % 
% % Input Variables
% % 
% % type is the boolean which selects A or C weighting.  
% %      0 chooses A-weighting.  
% %      1 chooses C-weighting.  default is 0, A-weighting
% % 
% % y is the multichannel input time record in (Pa).  Processsing assumes 
% %      that y has more channels than time record samples. 
% %      y=randn(10000, 10) is the default.  
% % 
% % Fs is the sampling rate in Hz.  default is 50000 Hz.
% %
% % settling_time is the time it takes the filter to settle (seconds).
% %      default is settling_time=0.1;
% % 
% % ***********************************************************
% % 
% % Output Variables
% % 
% % yAC is the A or C-weighted time record in (Pa).
% % 
% % errors indicates whether the resampled sampling rate is within the 
% %          tolerance of 5000 Hz.  Within tolerance 0 is output.  
% %          outside of tolerance 1 is output.   
% % 
% % ***********************************************************
%  
%  
% Example='1';
%  
% type=0;           % 0 selects A-weighting
%  
% y=randn(1, 50000);% (Pa) waveform
%                   % y should have the size (num_channels, num_datasample)
%                   
% Fs=50000;         % (Hz) Sampling rate 
%  
% settling_time=0.1;% (seconds) Time it takes the filter to settle.
% 
% [yA]=ACweight_time_filter(type, y, Fs, settling_time);
% 
% % Plot the results!
% t=1/Fs*(0:(Fs-1));
% figure(1); plot(t, y, 'k'); hold on; plot(t, yA, 'g'); legend('Linear Time Record', 'A-weighted');
% 
% type=1;           % 1 selects C-weighting
% [yC]=ACweight_time_filter(type, y, Fs, settling_time);
% 
% % Plot the results!
% figure(2); plot(t, y, 'k'); hold on; plot(t, yC, 'g'); legend('Linear Time Record', 'C-weighted');
% 
% % **********************************************************************
% % 
% % References
% % 
% % IEC/CD 1672: Electroacoustics-Sound Level Meters, Nov. 1996. 
% % 
% % ANSI S1.4: Specifications for Sound Level Meters, 1983. 
% % 
% % ***********************************************************
% % 
% % 
% % Subprograms
% %
% % This program requires the Matlab Signal Processing Toolbox
% %
% % 
% % List of Dependent Subprograms for 
% % ACweight_time_filter
% % 
% % 
% % Program Name   Author   FEX ID#
% % 1) ACdsgn		
% % 2) convert_double		
% % 3) estimatenoise		John D'Errico		16683	
% % 4) filter_settling_data		
% % 5) get_p_q		
% % 6) moving		Aslak Grinsted		8251	
% % 7) sub_mean		
% % 8) wsmooth		Damien Garcia		NA			
% % 
% % ***********************************************************
% % 
% % Written by Edward L. Zechmann   
% % 
% %     date    June        2007
% % 
% % modified 15 November    2007    Updated comments
% % 
% % modified 19 December    2007    Added a resampling routine 
% %                                 to improve accruacy at low and
% %                                 high frequencies 
% %                                 Updated comments
% %                                 
% % modified 15 August      2008    Modified resampling routine to use 
% %                                 a butterworth filter.  
% %                                 Updated comments
% % 
% % modified  9 September   2008    Reverted back to resample Matlab
% %                                 Program.  
% % 
% % modified 18 September   2008    Updated Comments
% % 
% % modified 18 September   2008    Modified resampling method
% % 
% % modified  6 December    2008    Added filter settling 
% % 
% % modified 10 December    2008    Added get_p_q to generalize the
% %                                 selection of the p and q for setting
% %                                 the resampling rate.  
% %  
% %                                 Combined Aweight_time_filter and
% %                                 Cweight_time_filter.  
% %  
% % modified 16 December    2008    Used convolution to make filter
% %                                 coefficients (b and a) into  
% %                                 arrays from cell arrays.
% %  
% % modified  5 January     2008    Added sub_mean to remove running
% %                                 average using a time constant set to 
% %                                 5 Hz frequency.
% %                                 
% % 
% % ***********************************************************
% % 
% % 
% % Please feel free to modify this code.
% % 
% % See also: adsgn, cdsgn, resample, filter_settling_data
% % 

if nargin < 1 || isempty(type) || ~isnumeric(type)
    type=0;
end

if nargin < 2 || isempty(y) || ~isnumeric(y)
    y=randn(10000, 10);
end

if nargin < 3 || isempty(Fs) || ~isnumeric(Fs)
    Fs=50000;
end

if (nargin < 4 || isempty(settling_time)) || ~isnumeric(settling_time)
    settling_time=0.1;
end

% set the flag for indicating whether to resample
flag0=0;


% Remove the running average from the signal.
% The time constant should be less than half the lowest frequency to
% resolve.
[y]=sub_mean(y, Fs, 5);


if Fs <= 120000
    p=1;
    q=1;
    Fsn=Fs;
    errors=0;
else
    tol=20000;
    [Fsn, p, q, errors]=get_p_q(Fs, 120000, tol);
end

if ~isequal(Fsn,Fs)
    flag0=1;
end

[m1 n1]=size(y);

% set the flag for indicating whether to transpose
flag1=0;

% Transpose the data so that the filters act along the time records 
% not along the channels.   
if m1 > n1
    flag1=1;
    y=y';
    [m1, n1]=size(y);
end

% Resample the data to a reasonable sampling rate to keep the A or 
% C-weighting filter as stable as possible.  
if isequal(flag0, 1)
    [y11, num_settle_pts]=filter_settling_data(Fs, y(1, :), settling_time);
    buf=resample([y11 y(1, :)], p, q);
    buf=buf(ceil((num_settle_pts*p/q)+1):end);
    
    n2=length(buf);
    y2=zeros(m1,n2);
    y2(1, :)=buf;
    clear('buf');
    if m1 > 1
        for e1=2:m1;
            [y11, num_settle_pts]=filter_settling_data(Fs, y(e1, :), settling_time);
            buf=resample([y11 y(e1, :)], p, q);
            y2(e1, :)=buf(ceil((num_settle_pts*p/q)+1):end);
        end
    end
else
    y2=y;
end


clear('y');

% Design the A or C-weighting filter

[Ba, Aa, Bc, Ac] = ACdsgn(Fsn);

if isequal(type, 0)
    B=Ba;
    A=Aa;
    %[B,A] = adsgn(Fsn);
else
    B=Bc;
    A=Ac;
    %[B,A] = cdsgn(Fsn);
end


yAC=zeros(size(y2));
   
% Apply the A or C-weighting filter
for e1=1:m1;
    [y11, num_settle_pts]=filter_settling_data(Fsn, y2(e1, :), settling_time);
    buf=[y11 y2(e1, :)];
    buf = real(filter(B, A, buf ));
    yAC(e1, :)=buf((num_settle_pts+1):end);
end




% Resample if necessary so the output has the same size as the input 
if isequal(flag0, 1)
    [y11, num_settle_pts]=filter_settling_data(Fs, yAC(1, :), settling_time);
    buf=resample([y11 yAC(1, :)], q, p);
    buf=buf(ceil((num_settle_pts*q/p)+1):end);
    
    n2=length(buf);
    yAACC=zeros(m1,n2);
    yAACC(1, :)=buf;
    clear('buf');
    if m1 > 1
        for e1=2:m1;
            [y11, num_settle_pts]=filter_settling_data(Fs, yAC(e1, :), settling_time);
            buf=resample([y11 yAC(e1, :)], q, p);
            yAACC(e1, :)=buf(ceil((num_settle_pts*q/p)+1):end);
        end
    end
else
    yAACC=yAC;
end


clear('yAC');

% Append the data to the output variable
if isequal(flag0, 1)
    if n2 > n1
        n2=n1;
    end

    yAC=yAACC(1:m1, 1:n2);
else
    yAC=yAACC(1:m1, 1:n1);
end

% Output yAC must have the same size as the original input matrix y.
% Transpose yAC if necessary. 
if isequal(flag1, 1)
    yAC=yAC';
end