WARPEDFILTERS - Frequency-warped band-limited filters

Usage

[g,a,fc]=warpedfilters(freqtoscale,scaletofreq,fs,fmin,fmax,bins,Ls);

Input parameters

Output parameters

Description

[g,a,fc]=warpedfilters(freqtoscale,scaletofreq,fs,fmin,fmax,bins,Ls) constructs a set of band-limited filters g which cover the required frequency range fmin-fmax with bins filters per scale unit. The filters are always centered at full (fractional \(k/bins\)) scale units, where the first filter is selected such that its center is lower than fmin.

By default, a Hann window on the frequency side is choosen, but the window can be changed by passing any of the window types from firwin as an optional parameter. Run getfield(getfield(arg_firwin,'flags'),'wintype') to get a cell array of window types available.

With respect to the selected scale, all filters have equal bandwidth and are uniformly spaced on the scale axis, e.g. if freqtoscale is \(\log(x)\), then we obtain constant-Q filters with geometric spacing. The remaining frequency intervals not covered by these filters are captured one or two additional filters (high-pass always, low-pass if necessary). The signal length Ls is required in order to obtain the optimal normalization factors.

Attention: When using this function, the user needs to be aware of a number of things:

1. Although the freqtoscale and scaletofreq can be chosen freely, it is assumed that freqtoscale is an invertible, increasing function from \(\mathbb{R}\) or \(\mathbb{R}^+\) onto \(\mathbb{R}\) and that freqtoscale is the inverse function.
2. If freqtoscale is from \(\mathbb{R}^+\) onto \(\mathbb{R}\), then necessarily \(freqtoscale(0) = -\infty\).
3. If the slope of freqtoscale is (locally) too steep, then there is the chance that some filters are effectively \(0\) or have extremely low bandwidth (1-3 samples), and consequently very poor localization in time. If freqtoscale is from \(\mathbb{R}^+\) onto \(\mathbb{R}\) then this usually occurs close to the DC component and can be alleviated by increasing fmin.
4. Since the input parameter bins is supposed to be integer, freqtoscale and scaletofreq have to be scaled appropriately. Note that \(freqtoscale(fs)\) is in some sense proportional to the resulting number of frequency bands and inversely proportional to the filter bandwidths. For example, the ERB scale defined by \(21.4\log_{10}(1+f/228.8)\) works nicely out of the box, while the similar mel scale \(2595\log_{10}(1+f/700)\) most likely has to be rescaled in order not to provide a filter bank with 1000s of channels.

If any of these guidelines are broken, this function is likely to break or give undesireable results.

By default, a Hann window is chosen as the transfer function prototype, but the window can be changed by passing any of the window types from firwin as an optional parameter.

The integer downsampling rates of the channels must all divide the signal length, filterbank will only work for input signal lengths being multiples of the least common multiple of the downsampling rates. See the help of filterbanklength. The fractional downsampling rates restrict the filterbank to a single length L=Ls.

[g,a]=warpedfilters(...,'regsampling') constructs a non-uniform filterbank with integer subsampling factors.

[g,a]=warpedfilters(...,'uniform') constructs a uniform filterbank where the the downsampling rate is the same for all the channels. This results in most redundant representation, which produces nice plots.

[g,a]=warpedfilters(...,'fractional') constructs a filterbank with fractional downsampling rates a. This results in the least redundant system.

[g,a]=warpedfilters(...,'fractionaluniform') constructs a filterbank with fractional downsampling rates a, which are uniform for all filters except the "filling" low-pass and high-pass filters can have different fractional downsampling rates. This is usefull when uniform subsampling and low redundancy at the same time are desirable.

The filters are intended to work with signals with a sampling rate of fs.

warpedfilters accepts the following optional parameters:

'bwmul',bwmul

Bandwidth variation factor. Multiplies the calculated bandwidth. Default value is 1. If the value is less than one, the system may no longer be painless.

'complex'

Construct a filterbank that covers the entire frequency range. When missing, only positive frequencies are covered.

'redmul',redmul

Redundancy multiplier. Increasing the value of this will make the system more redundant by lowering the channel downsampling rates. Default value is 1. If the value is less than one, the system may no longer be painless.

Examples:

In the first example, we use the ERB scale functions freqtoerb and erbtofreq to construct a filter bank and visualize the result:

[s,fs] = gspi; % Get a test signal
Ls = numel(gspi);

% Fix some parameters
fmax = fs/2;
bins = 1;

% Compute filters, using fractional downsampling
[g,a,fc]=warpedfilters(@freqtoerb,@erbtofreq,fs,0,fmax,bins,...
                       Ls,'bwmul',1.5,'real','fractional');

% Plot the filter transfer functions
figure(1);
filterbankfreqz(g,a,Ls,'plot','linabs','posfreq');
title('ERBlet filter transfer functions');

% Compute the frame bounds
gf=filterbankresponse(g,a,Ls,'real'); framebound_ratio = max(gf)/min(gf);
disp(['Painless system frame bound ratio of ERBlets: ',...
     num2str(framebound_ratio)]);

% Plot the filter bank coefficients of the test signal
figure(2);
c=filterbank(s,g,a);
plotfilterbank(c,a,fc,fs,60);
title('ERBlet transform of the test signal');

In the second example, we look at the same test signal using a constant-Q filter bank with 4 bins per scale unit and the standard (semi-regular) sampling scheme:

[s,fs] = gspi; % Get a test signal
Ls = numel(gspi);

% Fix some parameters
fmax = fs/2;
bins = 1;

% Define the frequency-to-scale and scale-to-frequency functions
warpfun_log = @(x) 10*log(x);
invfun_log = @(x) exp(x/10);

bins_hi = 4; % Select bins/unit parameter
fmin = 50; % The logarithm's derivative 1/x tends to Inf for x towards 0

% Compute filters, using fractional downsampling
[g,a,fc]=warpedfilters(warpfun_log,invfun_log,fs,fmin,fmax,bins_hi,Ls,'bwmul',1,'real');

% Plot the filter transfer functions
figure(1);
filterbankfreqz(g,a,Ls,'plot','linabs','posfreq');
title('constant-Q filter transfer functions (4 bins)');

% Compute the frame bounds
gf=filterbankresponse(g,a,Ls,'real'); framebound_ratio = max(gf)/min(gf);
disp(['Painless system frame bound ratio (constant-Q - 4 bins): ', num2str(framebound_ratio)]);

% Plot the filter bank coefficients of the test signal
figure(2);
c=filterbank(s,g,a);
plotfilterbank(c,a,fc,fs,60);
title('constant-Q transform of the test signal (4 bins)');