matlab中MakeERBFilters源代码

以下是MATLAB中的MakeERBFilters源代码：

function [filters, freqs] = MakeERBFilters(fs, numChannels, lowFreq)
%[filters, freqs] = MakeERBFilters(fs, numChannels, lowFreq)
%
%Create a set of Gammatone filters spaced according to the Equivalent
%Rectangular Bandwidth ERB. Filters are spaced according to the equation:
%
%   ERB = cF
%
%where F is the frequency and c is a constant. Here, c=1.0. Each filter is
%defined by the difference equation:
%
%   y[n] = b0*h0*x[n] - a1*y[n-1] - a2*y[n-2]
%
%where x is the input signal, y is the output signal, and h0, b0, a1, and a2
%are constants that define the filter.
%
%ARGUMENTS:
%   fs          -   the sampling rate in Hz.
%   numChannels -   the number of filters to create.
%   lowFreq     -   the frequency of the lowest filter in Hz.
%
%OUTPUTS:
%   filters     -   a matrix of shape (numChannels, filterLength) where
%                   filterLength is the number of samples in each filter.
%   freqs       -   a vector of shape (numChannels,) containing the center
%                   frequency of each filter.
%
%Author: James Lyons 2012

T = 1/fs;
% Change the following parameters if you wish to use a different ERB scale
% Minimum ERB (Hz):
minBW = 25;
% Maximum ERB (Hz):
maxBW = fs/2;
% Set up the centre frequencies.
cf_array=zeros(1,numChannels);
% Centre frequencies based on the hoerl approximation to the ERB formula.
cf=(erbspace(minBW/4,maxBW,numChannels));
cf_array(:)=cf;
% Now create a set of filters to cover the required bandwidth.
% Each loop creates a triangular filter with a bandwidth
% defined by ERB, then convolves it with the previous
% filter to create the next. First filter is a lowpass
% filter to ensure we don't alias when we downsample.
filters=zeros(numChannels,round(4*T*fs));
filterLength = size(filters,2);
% Compute constant coefficients for gammatone filter:
b=zeros(numChannels,3);
a=zeros(numChannels,2);
gain=zeros(numChannels,1);
for i=1:numChannels,
    cf = cf_array(i);
    ERB = 24.7*(4.37*cf/1000+1); % 24.7 is a constant scaling factor
    B = 1.019*2*pi*ERB; % 1.019 is the bandwidth in ERB.
    r = exp(-B*T);
    theta = 2*pi*cf*T;
    b(i,1) = T*T*2*pi*cf/B;
    b(i,2) = 0;
    b(i,3) = -T*T*2*pi*cf/B;
    a(i,1) = -2*r*cos(theta);
    a(i,2) = r*r;
    gain(i) = abs( (exp(1i*theta)-r)*T ).^2 / (2*pi*cf);
end
% First filter is lowpass:
filters(1,:) = gammatone_lowpass(b(1,:), a(1,:), gain(1), filterLength);
% Remaining filters are created by convolving previous filter with
% triangle filter:
for i=2:numChannels,
    ERB = 24.7*(4.37*cf_array(i)/1000+1); % 24.7 is a constant scaling factor
    B = 1.019*2*pi*ERB; % 1.019 is the bandwidth in ERB.
    filters(i,:) = gammatone_triangle(b(i,:), a(i,:), gain(i), B, filters(i-1,:));
end
% Finally, apply the Hanning window to each filter:
win = hanning(filterLength)';
filters = filters.*repmat(win,numChannels,1);
% compute the frequency response of each filter for plotting purposes:
freqs = cf_array;

该函数用于创建一组按等效矩形带宽（ERB）间隔的Gammatone滤波器。滤波器由差分方程定义，并且根据其中心频率，每个滤波器的带宽在等效矩形带宽（ERB）尺度下等间距。函数的输入包括采样率、要创建的滤波器数量以及最低滤波器的频率。函数的输出包括每个滤波器的系数和频率响应。