nth_oct_time_filter2.m

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

function [fc_out, SP_levels, SP_peak_levels, SP_bands]=Nth_oct_time_filter2(SP, Fs, num_x_filter, N, fc, sensor, settling_time, filter_program)
% % Nth_oct_time_filter2: Calculates the Nth octave center frequencies, sound levels, peak levels, and time records
% % 
% % Syntax:
% % 
% % [fc_out, SP_levels, SP_peak_levels, SP_bands]=Nth_oct_time_filter2(SP, Fs, num_x_filter, N, fc, sensor, settling_time, filter_program);
% %
% % **********************************************************************
% % 
% % Description
% %
% % This program applies Nth octave band filters to the input time record.  
% % The program outputs the center frequency bands, the time average rms 
% % values, the peak values, and band filtered time records for each 
% % Nth octave band respectively.  
% % 
% % To optimize filter stability, the resample program
% % is used iteratively to make the sampling rate reasonable
% % before applying the third octave filters.  Using the resmaple program
% % iteratively improves the antialiasing ability of the resmple program.  
% % 
% % Nth_octdsgn computes the filter coefficients using a 3rd order
% % butterworth filter for an Nth octave band filter according to 
% % ANSI S1.11.     
% % 
% % To avoid phase shift, the filtfilt Matlab program can
% % be used to implement the one-Nth octave filters.  
% % 
% % The input and output variables are described in more detail in the
% % sections below respectively.  
% % 
% % **********************************************************************
% % 
% % Input Variables
% %
% % SP=randn(10, 50000);
% %                     % (Pa) is the time record of the sound pressure
% %                     % default is SP=rand(1, 50000);
% %
% % Fs=50000;           % (Hz) is the sampling rate of the time record.
% %                     % default is Fs=50000; Hz.
% %
% % num_x_filter=2;     % This is the number of times the time record 
% %                     % should be filtered.  
% %                     % default is num_x_filter=2;
% % 
% % N=3;                % is the number of frequency bands per octave.  
% %                     % Can be any number > 0.  
% %                     % Default is 3 for third octave bands.  
% % 
% % fc=[200, 250];      % (Hz) is an array of center frequencies for the third
% %                     % octave band filters
% %                     %
% %                     % if empty the third octave bands from 20 to 20000 Hz
% %                     % are used
% %                     %
% %                     % default is fc = [  20,  25, 31.5, 40,  50,  63,  ...
% %                     %  80,  100, 125,  160, 200, 250, 315, 400, 500, ...
% %                     %  630, 800, 1000, 1250, 1600, 2000, 2500, 3150, ...
% %                     %  4000, 5000, 6300, 8000, 10000, ...
% %                     %  12500, 16000, 20000];
% %
% % sensor=1;           % Constant integer input for selecting the sensor type
% %                     % 1 is for acoustic microphone Pref=20E-6 (Pa)
% %                     %
% %                     % 2 is for accelerometer output is in same
% %                     %   units as the input (m/s^2)
% %                     %
% %                     % 3 generic sensor multiply by 1: output is in same
% %                     %   units as the input
% %                     %
% %                     % default is sensor=1; For a microphone
% % 
% % settling_time=0.1;  % (seconds) Time requiered for the filter to settle 
% %                     % usually 0.1 seconds or less.  
% %                     % This quantity is usually frequency dependent.  
% %
% % filter_program=1;   % 1 is for using the filter progam otherwise the
% %                     % filtfilt program is used.  
% %                     % filter.m runs faster and may settle 
% %                     % more quickly.    
% %                     % filtfilt.m is used to remove phase shift. 
% %                     % default is filter_program=1 using filter progam.
% % 
% % **********************************************************************
% %
% % Output Variables
% %
% % fc_out          % (Hz) array of center frequencies
% %
% % SP_levels       % (units) sound pressure or other sensor metric for 
% %                 % each channel and each frequency band
% % 
% % SP_peak_levels  % (units) Maximum of the absolute value of the peak
% %                 % values for each frequency band
% %
% % SP_bands        % Time record for each mic channel and for each
% %                 % frequency band after filtering
% %
% %
% % **********************************************************************
% %
%
% Example='1';
%
% SP=randn(1, 50000);   % SP is the data variables in linear units such as
%                       % (Pa)
%
% Fs=50000;             % (Hz) Sampling rate
%
% num_x_filter=2;       % Number of times to filter the data.  Minimum 
%                       % value is 1.  Typically a value of 2 to 10 at low
%                       % frequencies (Fc < 100), num_x_filter=10 has a
%                       % significant phase shift.
%
% N=3;                  % Number of bands per octave.  
%
% fc=[];                % (Hz) Center frequency of the third-octave band
%
% sensor=1;             % acoustic microphone
%                       % output is in dB
% 
% settling_time=1;      % (seconds) Time requiered for the filter to settle 
%                       % usually 0.1 seconds or less.  
%                       % This quantity is usually frequency dependent.  
%
% filter_program=1;     % 1 is for using the filter progam otherwise the
%                       % filtfilt program is used.  
%                       % default is filter_program=1 using filter progam.
%
% [fc_out, SP_levels, SP_peak_levels, SP_bands]=Nth_oct_time_filter2(SP, Fs, num_x_filter, N, fc, sensor, settling_time, filter_program);
% 
% %
%
% Example='1';
% 
% % Compare the spectra of white noise, pink noise, and brown noise.  
% % 
% 
% x1 = spatialPattern([1,500000],0);    % white noise has a linearly
%                                       % increasing spectrum
%
% x2 = spatialPattern([1,500000],-1);   % pink noise has a constant
%                                       % spectrum
% 
% x3 = spatialPattern([1,500000],-2);   % brown noise has a linearly
%                                       % increasing spectra
%
% Fs=50000;         % (Hz) Sampling rate
%
% num_x_filter=2;   % Number of times to filter the data.  Minimum value is 1
%                   % typically a value of 2 to 10 at low
%                   % frequencies (Fc < 100), num_x_filter=10 has a
%                   % significant phase shift when using filter.
% 
% N=3;              % number of bands per octave.  
% 
% min_f=20;         % is the minimum frequency band to calculate (Hz).   
% max_f=20000;      % is the maximum frequency band to calculate (Hz).  
% 
% [fc] = nth_freq_band(N, min_f, max_f);
%
% sensor=1;         % acoustic microphone
%                   % output is in dB
% 
% settling_time=1;  % (seconds) Time requiered for the filter to settle 
%                   % usually 0.1 seconds or less.  
%                   % This quantity is usually frequency dependent.  
%
% filter_program=1; % 1 is for using the filter progam otherwise the
%                   % filtfilt program is used.  
%                   % default is filter_program=1 using filter progam.
%
% [fc_out1, SP_levels1]=Nth_oct_time_filter(x1, Fs, num_x_filter, N, min_f, max_f, sensor, settling_time, filter_program);
% [fc_out2, SP_levels2]=Nth_oct_time_filter(x2, Fs, num_x_filter, N, min_f, max_f, sensor, settling_time, filter_program);
% [fc_out3, SP_levels3]=Nth_oct_time_filter(x3, Fs, num_x_filter, N, min_f, max_f, sensor, settling_time, filter_program);
%
% % Plot the results
% figure(1); 
% semilogx(fc_out1, SP_levels1, 'color', [1 1 1],         'linewidth', 2,                    'marker', 's', 'MarkerSize', 8);
% hold on;
% semilogx(fc_out2, SP_levels2, 'color', [1 0.6 0.784],   'linewidth', 2, 'linestyle', '--', 'marker', 'o', 'MarkerSize', 8);
% semilogx(fc_out3, SP_levels3, 'color', [0.682 0.467 0], 'linewidth', 2, 'linestyle', ':',  'marker', 'x', 'MarkerSize', 12);
% set(gca, 'color', 0.7*[1 1 1]);
% legend({'White Noise', 'Pink Noise', 'Brown Noise'}, 'location', 'SouthEast');
% xlabel('Frequency Hz', 'Fontsize', 28);
% ylabel('Sound Pressure Level (dB ref. 20 \mu Pa)', 'Fontsize', 28);
% title('Classical Third Octave Band Spectra', 'Fontsize', 40);
% set(gca, 'Fontsize', 20);
% 
% % **********************************************************************
% % 
% % References
% % 
% % 1)  ANSI S1.11-1986 American National Stadard Specification for 
% %                     Octave-Band and Fractional-Octave-Band Analog 
% %                     and Digital Filters.
% % 
% % 
% % **********************************************************************
% % 
% % Subprograms
% %
% % This program requires the Matlab Signal Processing Toolbox
% % This program uses a recreation of  oct3dsgn	by Christophe Couvreur	69
% % 
% % 
% % 
% % List of Dependent Subprograms for 
% % Nth_oct_time_filter2
% % 
% % 
% % Program Name   Author   FEX ID#
% % 1) convert_double		
% % 2) estimatenoise		John D'Errico		16683	
% % 3) filter_settling_data		
% % 4) moving		Aslak Grinsted		8251	
% % 5) Nth_octdsgn		
% % 6) sub_mean		
% % 7) wsmooth		Damien Garcia							
% %
% % **********************************************************************
% % 
% % Program was written by Edward L. Zechmann
% % 
% %     date  7 December    2008    Copied content of Nth_oct_time_filter
% %                                 and modified code to use fc as an input
% %                                 variable.
% % 
% % modified  8 December    2008    Updated Comments
% % 
% % modified 10 December    2008    Updated Comments
% % 
% % modified 16 December    2008    Generlaized program to Nth Octave Bands
% % 
% % modified 22 December    2008    Updated Comments.  Finished Upgrade
% % 
% % modified  5 January     2008    Added sub_mean to remove running
% %                                 average using a time constant at one- 
% %                                 half the lowest center frequency.
% %                                 
% % 
% % **********************************************************************
% % 
% % Please feel free to modify this code.
% %
% % See Also: Nth_oct_time_filter, octave, resample, filter, filtfilt
% %

if (nargin < 1 || isempty(SP)) || ~isnumeric(SP)
    SP=rand(1, 50000);
end

% Make the data have the correct data type and size
[SP]=convert_double(SP);

[num_mics, num_pts]=size(SP);

if num_mics > num_pts
    SP=SP';
    [num_mics num_pts]=size(SP);
end

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

if (nargin < 3 || isempty(num_x_filter)) || ~isnumeric(num_x_filter)
    num_x_filter=2;
end

if (nargin < 4 || (isempty(N)) || ~isnumeric(N))
    N=20;
end

if (nargin < 5 || isempty(fc)) || ~isnumeric(fc)
    fc = [  20,  25, 31.5, 40,  50,  63,  80,  100,  125,  160, ...
        200, 250, 315, 400, 500, 630, 800, 1000, 1250, 1600, ...
        2000, 2500, 3150, 4000, 5000, 6300, 8000, 10000, ...
        12500, 16000, 20000];
end

if (nargin < 6 || isempty(sensor)) || ~isnumeric(sensor)
    sensor=1;
end

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

if (nargin < 8 || (isempty(filter_program)) || ~isnumeric(filter_program))
    filter_program=1;
end

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


% By hard coding, a 3rd order butterworth filter is used.  
n=3;

% Make sure center frequencies are positive.
fc=fc(fc > 0);

% sort the center frequencies to categorize teh center freuncies into
% ranges for resampling.
[fc ix]=sort(fc);
fc_out=fc;

% count the frequency bands
num_bands=length(fc);

% srr is the sampling rate ratio
srr=Fs./fc;

% initialize the regime to ones.
regime=zeros(num_bands, 1);

% Using a small range size keeps 
% increases the likelihood that the filters will be stable.  
% In testing, using range_size=2, 4, or 8 resulted in nearly identical
% sound pressure levels, however; range_size=2, is more liekly to result in 
% a stable filter.  
range_size=2;

% 
% For range_size=4; the regime follows the pattern
% 
%   srr range        Ouput type        rf value   regime number
%       srr <=    2  output zeros                      1
%   2 < srr <=    5  upsample          rf=4            2
%   5 < srr <=   20  do nothing        rf=1            3
%  20 < srr <=   80  downsample        rf=1/4          4
%  80 < srr <=  320  downsample        rf=1/16         5
% 320 < srr <= 1280  downsample        rf=1/64         6
%                      ...
%                      ...  iterative downsampling until the frequency
%                      ...  resolution limit is exceeded
%                      ...
%  fc <= 4/T         output zeros                      0