dicklyon@465: function [complex_transfns_freqs, ...
dicklyon@465:   stage_numerators, stage_denominators] = ...
dicklyon@465:   CARFAC_Transfer_Functions(CF, freqs, to_channels, from_channels)
dicklyon@465: % function [complex_transfns_freqs, ...
dicklyon@465: %   stage_z_numerators, stage_z_denominators] = ...
dicklyon@465: %   CARFAC_Transfer_Functions(CF, freqs, to_channels, from_channels)
dicklyon@465: % Return transfer functions as polynomials in z (nums & denoms);
dicklyon@465: % And evaluate them at freqs if it's given, to selected output,
dicklyon@465: %   optionally from selected starting points (from 0, input, by default).
dicklyon@465: %   complex_transfns_freqs has a row of complex gains per to_channel.
dicklyon@465: 
dicklyon@465: % always start with the rational functions, whether we want to return
dicklyon@465: % them or not:
dicklyon@465: [stage_numerators, stage_denominators] = CARFAC_Rational_Functions(CF);
dicklyon@465: 
dicklyon@465: if nargin >= 2
dicklyon@465:   % Evaluate at the provided list of frequencies.
dicklyon@465:   if ~isrow(freqs)
dicklyon@465:     if iscolumn(freqs)
dicklyon@465:       freqs = freqs';
dicklyon@465:     else
dicklyon@465:       error('Bad freqs_row in CARFAC_Transfer_Functions');
dicklyon@465:     end
dicklyon@465:   end
dicklyon@465:   if any(freqs < 0)
dicklyon@465:     error('Negatives in freqs_row in CARFAC_Transfer_Functions');
dicklyon@465:   end
dicklyon@465:   z_row = exp((i * 2 * pi / CF.fs) * freqs);  % z = exp(sT)
dicklyon@465:   gains = Rational_Eval(stage_numerators, stage_denominators, z_row);
dicklyon@465:   
dicklyon@465:   % Now multiply gains from input to output places; use logs?
dicklyon@465:   log_gains = log(gains);
dicklyon@465:   cum_log_gains = cumsum(log_gains);
dicklyon@465:   
dicklyon@465:   % And figure out which cascade products we want:
dicklyon@465:   n_ch = CF.n_ch;
dicklyon@465:   if nargin < 3
dicklyon@465:     to_channels = 1:n_ch;
dicklyon@465:   end
dicklyon@465:   if isempty(to_channels) || any(to_channels < 1 | to_channels > n_ch)
dicklyon@465:     error('Bad to_channels in CARFAC_Transfer_Functions');
dicklyon@465:   end
dicklyon@465:   if nargin < 4 || isempty(from_channels)
dicklyon@465:     from_channels = 0;  % tranfuns from input, called channel 0.
dicklyon@465:   end
dicklyon@465:   if length(from_channels) == 1
dicklyon@465:     from_channels = from_channels * ones(length(to_channels));
dicklyon@465:   end
dicklyon@465:   % Default to cum gain of 1 (log is 0), from input channel 0:
dicklyon@465:   from_cum = zeros(length(to_channels), length(z_row));
dicklyon@465:   not_input = from_channels > 0;
dicklyon@465:   from_cum(not_input, :) = cum_log_gains(from_channels(not_input), :);
dicklyon@465:   log_transfns = cum_log_gains(to_channels, :) - from_cum;
dicklyon@465:   complex_transfns_freqs = exp(log_transfns);
dicklyon@465: else
dicklyon@465:   % If no freqs are provided, do nothing but return the stage info above:
dicklyon@465:   complex_transfns_freqs = [];
dicklyon@465: end
dicklyon@465: 
dicklyon@465: 
dicklyon@465: 
dicklyon@465: function gains = Rational_Eval(numerators, denominators, z_row)
dicklyon@465: % function gains = Rational_Eval(numerators, denominators, z_row)
dicklyon@465: % Evaluate rational function at row of z values.
dicklyon@465: 
dicklyon@465: zz = [z_row .* z_row; z_row; ones(size(z_row))];
dicklyon@465: % dot product of each poly row with each [z2; z; 1] col:
dicklyon@465: gains = (numerators * zz) ./ (denominators * zz);
dicklyon@465: 
dicklyon@465: 
dicklyon@465: function [stage_numerators, stage_denominators] = ...
dicklyon@465:   CARFAC_Rational_Functions(CF, chans)
dicklyon@465: % function [stage_z_numerators, stage_z_denominators] = ...
dicklyon@465: %   CARFAC_Rational_Functions(CF, chans)
dicklyon@465: % Return transfer functions of all stages as rational functions.
dicklyon@465: 
dicklyon@465: if nargin < 2
dicklyon@465:   n_ch = CF.n_ch;
dicklyon@465:   chans = 1:n_ch;
dicklyon@465: else
dicklyon@465:   n_ch = length(chans);
dicklyon@465: end
dicklyon@465: coeffs = CF.filter_coeffs;
dicklyon@465: r = coeffs.r_coeffs(chans);
dicklyon@465: a = coeffs.a_coeffs(chans) .* r;
dicklyon@465: c = coeffs.c_coeffs(chans) .* r;
dicklyon@465: r2 = r .* r;
dicklyon@465: h = coeffs.h_coeffs(chans);
dicklyon@465: g = coeffs.g_coeffs(chans);
dicklyon@465: stage_denominators = [ones(n_ch, 1), -2 * a, r2];
dicklyon@465: stage_numerators = [g .* ones(n_ch, 1), g .* (-2 * a + h .* c), g .* r2];
dicklyon@465: 
dicklyon@465: 
dicklyon@465: %% example
dicklyon@465: % CF = CARFAC_Design
dicklyon@465: % f = (0:100).^2;  % frequencies to 10 kHz, unequally spaced
dicklyon@465: % to_ch = 10:10:96;  % selected output channels
dicklyon@465: % from_ch = to_ch - 10;  % test the inclusion of 0 in from list
dicklyon@465: % tf = CARFAC_Transfer_Functions(CF, f, to_ch, from_ch);
dicklyon@465: % figure
dicklyon@465: % plot(f, 20*log(abs(tf)')/log(10));  % dB vs lin. freq for 10 taps
dicklyon@465: