dicklyon@527: % Copyright 2012, Google, Inc. dicklyon@527: % Author: Richard F. Lyon dicklyon@527: % dicklyon@527: % This Matlab file is part of an implementation of Lyon's cochlear model: dicklyon@527: % "Cascade of Asymmetric Resonators with Fast-Acting Compression" dicklyon@527: % to supplement Lyon's upcoming book "Human and Machine Hearing" dicklyon@527: % dicklyon@527: % Licensed under the Apache License, Version 2.0 (the "License"); dicklyon@527: % you may not use this file except in compliance with the License. dicklyon@527: % You may obtain a copy of the License at dicklyon@527: % dicklyon@527: % http://www.apache.org/licenses/LICENSE-2.0 dicklyon@527: % dicklyon@527: % Unless required by applicable law or agreed to in writing, software dicklyon@527: % distributed under the License is distributed on an "AS IS" BASIS, dicklyon@527: % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. dicklyon@527: % See the License for the specific language governing permissions and dicklyon@527: % limitations under the License. dicklyon@527: dicklyon@526: function [complex_transfns_freqs, ... dicklyon@533: stage_numerators, stage_denominators, group_delays] = ... dicklyon@533: CARFAC_Transfer_Functions(CF, freqs, to_channels, from_channels) dicklyon@526: % function [complex_transfns_freqs, ... dicklyon@533: % stage_numerators, stage_denominators, group_delays] = ... dicklyon@533: % CARFAC_Transfer_Functions(CF, freqs, to_channels, from_channels) dicklyon@533: % dicklyon@526: % Return transfer functions as polynomials in z (nums & denoms); dicklyon@526: % And evaluate them at freqs if it's given, to selected output, dicklyon@526: % optionally from selected starting points (from 0, input, by default). dicklyon@526: % complex_transfns_freqs has a row of complex gains per to_channel. dicklyon@526: dicklyon@526: % always start with the rational functions, whether we want to return dicklyon@526: % them or not: dicklyon@526: [stage_numerators, stage_denominators] = CARFAC_Rational_Functions(CF); dicklyon@526: dicklyon@526: if nargin >= 2 dicklyon@526: % Evaluate at the provided list of frequencies. dicklyon@526: if ~isrow(freqs) dicklyon@526: if iscolumn(freqs) dicklyon@526: freqs = freqs'; dicklyon@526: else dicklyon@526: error('Bad freqs_row in CARFAC_Transfer_Functions'); dicklyon@526: end dicklyon@526: end dicklyon@526: if any(freqs < 0) dicklyon@526: error('Negatives in freqs_row in CARFAC_Transfer_Functions'); dicklyon@526: end dicklyon@526: z_row = exp((i * 2 * pi / CF.fs) * freqs); % z = exp(sT) dicklyon@526: gains = Rational_Eval(stage_numerators, stage_denominators, z_row); dicklyon@526: dicklyon@526: % Now multiply gains from input to output places; use logs? dicklyon@526: log_gains = log(gains); dicklyon@533: cum_log_gains = cumsum(log_gains); % accum across cascaded stages dicklyon@526: dicklyon@526: % And figure out which cascade products we want: dicklyon@526: n_ch = CF.n_ch; dicklyon@526: if nargin < 3 dicklyon@526: to_channels = 1:n_ch; dicklyon@526: end dicklyon@526: if isempty(to_channels) || any(to_channels < 1 | to_channels > n_ch) dicklyon@526: error('Bad to_channels in CARFAC_Transfer_Functions'); dicklyon@526: end dicklyon@526: if nargin < 4 || isempty(from_channels) dicklyon@526: from_channels = 0; % tranfuns from input, called channel 0. dicklyon@526: end dicklyon@526: if length(from_channels) == 1 dicklyon@528: from_channels = from_channels * ones(1,length(to_channels)); dicklyon@526: end dicklyon@526: % Default to cum gain of 1 (log is 0), from input channel 0: dicklyon@526: from_cum = zeros(length(to_channels), length(z_row)); dicklyon@526: not_input = from_channels > 0; dicklyon@526: from_cum(not_input, :) = cum_log_gains(from_channels(not_input), :); dicklyon@526: log_transfns = cum_log_gains(to_channels, :) - from_cum; dicklyon@526: complex_transfns_freqs = exp(log_transfns); dicklyon@533: dicklyon@533: if nargout >= 4 dicklyon@533: phases = imag(log_gains); % no wrapping problem on single stages dicklyon@533: cum_phases = cumsum(phases); % so no wrapping here either dicklyon@533: group_delays = -diff(cum_phases')'; % diff across frequencies dicklyon@533: group_delays = group_delays ./ (2*pi*repmat(diff(freqs), n_ch, 1)); dicklyon@533: end dicklyon@526: else dicklyon@526: % If no freqs are provided, do nothing but return the stage info above: dicklyon@526: complex_transfns_freqs = []; dicklyon@526: end dicklyon@526: dicklyon@526: dicklyon@526: dicklyon@526: function gains = Rational_Eval(numerators, denominators, z_row) dicklyon@526: % function gains = Rational_Eval(numerators, denominators, z_row) dicklyon@526: % Evaluate rational function at row of z values. dicklyon@526: dicklyon@526: zz = [z_row .* z_row; z_row; ones(size(z_row))]; dicklyon@526: % dot product of each poly row with each [z2; z; 1] col: dicklyon@526: gains = (numerators * zz) ./ (denominators * zz); dicklyon@526: dicklyon@526: dicklyon@526: function [stage_numerators, stage_denominators] = ... dicklyon@528: CARFAC_Rational_Functions(CF) dicklyon@526: % function [stage_z_numerators, stage_z_denominators] = ... dicklyon@526: % CARFAC_Rational_Functions(CF, chans) dicklyon@526: % Return transfer functions of all stages as rational functions. dicklyon@526: dicklyon@528: n_ch = CF.n_ch; dicklyon@528: coeffs = CF.filter_coeffs; dicklyon@528: min_zeta = CF.filter_params.min_zeta; dicklyon@528: dicklyon@530: a0 = coeffs.a0_coeffs; dicklyon@530: c0 = coeffs.c0_coeffs; dicklyon@530: zr = coeffs.zr_coeffs; dicklyon@528: dicklyon@528: % get r, adapted if we have state: dicklyon@530: r = coeffs.r1_coeffs; dicklyon@528: if isfield(CF, 'filter_state') dicklyon@528: state = CF.filter_state; dicklyon@528: zB = state.zB_memory; % current extra damping dicklyon@530: r = r - zr .* zB; dicklyon@526: else dicklyon@528: zB = 0; dicklyon@526: end dicklyon@528: dicklyon@530: g = CARFAC_Stage_g(coeffs, zB); dicklyon@528: a = a0 .* r; dicklyon@528: c = c0 .* r; dicklyon@526: r2 = r .* r; dicklyon@528: h = coeffs.h_coeffs; dicklyon@528: dicklyon@526: stage_denominators = [ones(n_ch, 1), -2 * a, r2]; dicklyon@526: stage_numerators = [g .* ones(n_ch, 1), g .* (-2 * a + h .* c), g .* r2]; dicklyon@526: dicklyon@526: dicklyon@526: %% example dicklyon@526: % CF = CARFAC_Design dicklyon@526: % f = (0:100).^2; % frequencies to 10 kHz, unequally spaced dicklyon@526: % to_ch = 10:10:96; % selected output channels dicklyon@526: % from_ch = to_ch - 10; % test the inclusion of 0 in from list dicklyon@526: % tf = CARFAC_Transfer_Functions(CF, f, to_ch, from_ch); dicklyon@526: % figure dicklyon@526: % plot(f, 20*log(abs(tf)')/log(10)); % dB vs lin. freq for 10 taps dicklyon@526: