nth_oct_time_filter2.m

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

% Determine the regime for each frequency band
% 
for e1=1:num_bands;

    if srr(e1) <= 2
        % Output Zeros since Center Frequency exceeds the Nyquist Frequency
        regime(e1)=1;
    elseif  fc(e1) <= 4*Fs/(num_pts)
        % Output Zeros since Center Frequency is below the lower
        % frequency resolution limit
        regime(e1)=0;
    elseif (2 < srr(e1)) && logical(srr(e1)<= 5)
        % The signal is upsampled
        regime(e1)=2;
    else
        % The signal is downsampled
        regime(e1)=ceil(log(srr(e1)/5)./log(range_size))+2;   
    end

end



% set the reference sensor value
switch sensor

    case 1
        % reference sound pressure
        Pref=20*10^(-6); % Pa
        
    case 2
        % reference acceleration
        Pref=1; % m/s^2
        
    case 3
        Pref=1;
        
    otherwise
        Pref=1;
end


% Initialize the output variables to reasonable values
SP_levels=0;
SP_bands=0;
SP_peak_levels=0;
cregime=0;

% Initialize the previous channel index
pe1=0;



% Calculate the sound pressure levels, peak levels, and band time records
if logical(num_mics > 0) && ((logical(num_bands > 0) && logical(num_pts > 0)))


    SP_levels=zeros(num_mics, num_bands);
    SP_peak_levels=zeros(num_mics, num_bands);

    % The time record in each band is only initialized if
    % it is output.  This saves memory and reduces processing time.  
    if nargout > 2
        SP_bands=zeros(num_mics, num_bands, num_pts);
    end
    
    
    
    for e1=1:num_mics;

        for e2=num_bands:(-1):1;

            % Set the third octave center frequency
            fc2=fc(e2);

            % initialize the flag which determines whether to output zeros
            flag2=0;

            % Determine whether to output zeros.
            % Initialize the data variable SP_trunc2.
            if isequal(regime(e2), 0) || isequal(regime(e2), 1)
                flag2=1;
            else
                if  ~isequal(e1, pe1) || logical(cregime > regime(e2))
                    SP_trunc2=SP(e1, :);   
                else
                    if cregime < 3 && regime(e2) >=3
                        SP_trunc2=SP(e1, :); 
                    end
                end
                            
            end
            
            
            
            
            if isequal(flag2, 0)

                Fs2=Fs;
                
                % Resample the time record if necessary
                if regime(e2) > 1
                    
                    if isequal(regime(e2), 2)
                        
                        % set the resample rate
                        upsample_factor=5;
                        Fs2=Fs*upsample_factor;
                        
                        % upsample if necessary
                        if ~isequal(cregime, regime(e2)) || ~isequal(e1, pe1)
                            
                            % Determine appropriate data for settling
                            % the filter.
                            [y2, num_pts_se]=filter_settling_data(Fs, SP_trunc2, settling_time);
                            
                            % Resample the time record with the preappended
                            % settling data.  
                            SP_trunc2=resample([y2 SP_trunc2], upsample_factor, 1);
                            
                            % Remove the settling data from the time
                            % record. 
                            SP_trunc2=SP_trunc2((floor(num_pts_se*upsample_factor)+1):end);
                        end
                        
                    elseif regime(e2) > 3 
                        
                        Fs2=Fs/(range_size^(regime(e2)-3));
                        
                        % Down Sample if necessary
                        if ~isequal(cregime, regime(e2)) || ~isequal(e1, pe1)
                            
                            % determine the amount of down sampling
                            if cregime >= 3 && logical(regime(e2) > cregime)
                                num_iter=regime(e2)-cregime;
                            else
                                num_iter=regime(e2)-3;
                            end
                            
                            
                            
                            for e3=1:num_iter;
                                
                                % Calculate the sampling rate for selecting
                                % the filter settling data
                                Fs22=Fs/(range_size^(regime(e2)-3+e3-num_iter-1));
                                
                                % Determine appropriate data for settling
                                % the filter.
                                [y2, num_pts_se]=filter_settling_data(Fs22, SP_trunc2, settling_time);
                                
                                % Resample the time record with the preappended
                                % settling data.  
                                SP_trunc2=resample([y2 SP_trunc2], 1, range_size);
                                
                                % Remove the settling data from the time
                                % record.  
                                SP_trunc2=SP_trunc2((floor(num_pts_se/range_size)+1):end);
                                
                            end
                            
                        end

                    end
                    
                end
                
                % transfer the time record to a buffer variable
                SP_trunc22=SP_trunc2;

               
                % Calculate the 1/3 octave band filter coefficients
                [Bc, Ac]=Nth_octdsgn(Fs2, fc2, N, n);
                
                
                % Determine the size of data needed to keep the filter from failing  
                flag3=0;
                [num_chans, data_pts]=size(SP_trunc22);
                filtorder=max([size(Bc, 2), size(Ac, 2)]);
                
                % If there is too little data then the filter will fail
                % replicate the time record several times so that there is 
                % enough data to keep the filter from failing. 
                if data_pts < 5*filtorder
                    flag3=1;
                    num_reps=ceil(5*filtorder/data_pts);
                    SP_trunc22 = repmat(SP_trunc22, 1, num_reps);
                end
                
                % Determine appropriate data for settling
                % the filter.
                [y2, num_pts_se]=filter_settling_data(Fs2, SP_trunc22, settling_time);
                
                % Preappend the filter settling data.
                SP_trunc22=[y2 SP_trunc22];
                
                
                % Apply the 1/3 octave bandpass filter to the data.
                for e3=1:num_x_filter;
                    
                    if isequal(filter_program, 1)
                        SP_trunc22 = filter(Bc,Ac,SP_trunc22);
                    else
                        SP_trunc22 = filtfilt(Bc,Ac,SP_trunc22);
                    end
                    
                end
                
                
                % Remove the settling data from the time
                % record.  
                SP_trunc22=SP_trunc22((num_pts_se+1):end);
                
                % Remove any replicated data used to keep the filter from
                % failing.  
                if flag3 == 1
                    SP_trunc22=SP_trunc22(1:num_chans, 1:data_pts);
                end
                
                
                % Only concatenate the 1/nth octave band time records
                % if the output variable exists.
                if nargout > 2
                    
                    % resample the output time record to the original
                    % sampling rate
                    
                    if regime(e2) == 2 
                        % Determine appropriate data for settling
                        % the filter.
                        [y2, num_pts_se]=filter_settling_data(Fs*upsample_factor, SP_trunc22, settling_time);
                        
                        % Resample the time record to reduce the sampling
                        % rate.
                        SP_trunc22=resample([y2 SP_trunc22], 1, upsample_factor);
                        
                        % Remove the settling data from the time
                        % record.
                        SP_trunc22=SP_trunc22((floor(num_pts_se/upsample_factor)+1):end);
                        
                        
                    elseif regime(e2) > 3 
                        
                        % Calculate the down sampling factor applied to
                        % the data.
                        %
                        % This factor is called the reduction factor.  
                        reduction_factor=(range_size^(regime(e2)-3));
                        
                        % Determine appropriate data for settling
                        % the filter.
                        [y2, num_pts_se]=filter_settling_data(Fs/reduction_factor, SP_trunc22, settling_time);
                        
                        % Upsample the time record to the original 
                        % sampling rate.
                        SP_trunc22=resample([ y2 SP_trunc22], reduction_factor, 1);
                        
                        % Remove the settling data from the time
                        % record.
                        SP_trunc22=SP_trunc22((floor(num_pts_se*reduction_factor)+1):end);
                        
                    end
                    
                    % Append the time record for the microphone channel 
                    % to the output band. 
                    num_pts2=min([length(SP_trunc22) num_pts]);

                    for e3=1:num_pts2;
                        SP_bands(e1, ix(e2), e3)=SP_trunc22(e3);
                    end

                end
                
                
                % Calculate the levels and peak levels.  
                switch sensor

                    case 1
                        SP_levels(e1, ix(e2))=10*log10((norm(SP_trunc22)./sqrt(length(SP_trunc22))./Pref).^2);
                        SP_peak_levels(e1, ix(e2))=20*log10(max(abs(SP_trunc22))./Pref);
                    case 2
                        SP_levels(e1, ix(e2))=norm(SP_trunc22)./sqrt(length(SP_trunc22));
                        SP_peak_levels(e1, ix(e2))=max(abs(SP_trunc22));
                    case 3
                        SP_levels(e1, ix(e2))=norm(SP_trunc22)./sqrt(length(SP_trunc22));
                        SP_peak_levels(e1, ix(e2))=max(abs(SP_trunc22));
                    otherwise
                        SP_levels(e1, ix(e2))=norm(SP_trunc22)./sqrt(length(SP_trunc22));
                        SP_peak_levels(e1, ix(e2))=max(abs(SP_trunc22));
                end

            else
                
                %
                % If the center frequency is greater than the  
                % Nyquist freuqency or less than the minimum frequency 
                % limit then output a zero time record and indicate
                % that the levels and peak levels are null.
                % 

                if nargout > 3

                    for e3=1:num_pts2;
                        SP_bands(e1, ix(e2), e3)=0;
                    end

                end

                switch sensor

                    case 1
                        SP_levels(e1, ix(e2))=-10^15;
                        SP_peak_levels(e1, ix(e2))=-10^15;
                    case 2
                        SP_levels(e1, ix(e2))=0;
                        SP_peak_levels(e1, ix(e2))=0;
                    case 3
                        SP_levels(e1, ix(e2))=0;
                        SP_peak_levels(e1, ix(e2))=0;
                    otherwise
                        SP_levels(e1, ix(e2))=0;
                        SP_peak_levels(e1, ix(e2))=0;
                end
            end

            
            % Set the previous hcannel index.
            pe1=e1;
            
            % Set the previous regime number.
            cregime=regime(e2);
                
        end
    end
end