tomwalters@0: function [outputs, state1, state2] = PZBankStep2(inputs, pfreqs, pdamps, ...
tomwalters@0:     mindamp, maxdamp, xmin, xmax, rmin, rmax, za0, za1, za2, state1, state2, ...
tomwalters@0:     showconstellation, showtransfuns)
tomwalters@0: % function [outputs, state1, state2] = PZBankStep2(inputs, pfreqs, pdamps, ...
tomwalters@0: %     mindamp, maxdamp, xmin, xmax, rmin, rmax, za0, za1, za2, state1, state2, ...
tomwalters@0: %     showconstellation, showtransfuns)
tomwalters@0: % one time step of Pole-Zero-FilterBank, using row vectors of parameters
tomwalters@0: % inputs should contain shifted outputs from last time with new scalar
tomwalters@0: %   waveform input on the front:  inputs = [newin, outputs(1:N-1)];
tomwalters@0: % version 2 uses pole x and radius interpolation to adjust dampling
tomwalters@0: % the min and max are interpolation reference points, not limits to enforce
tomwalters@0: % and it should all work for multiple channels in parallel if the
tomwalters@0: % parameter vectors are all appropriately duplicated as well
tomwalters@0: 
tomwalters@0: if nargin < 15
tomwalters@0:     showconstellation = 0; % don't by default; value is used as fig num.
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if nargin < 16
tomwalters@0:     showtransfuns = 0; % don't by default; value is used as fig num.
tomwalters@0: end
tomwalters@0: 
tomwalters@0: [nstages, nchans] = size(pfreqs); % number of filter channels and stereo channels
tomwalters@0: 
tomwalters@0: damprate = 1/(maxdamp - mindamp); % to save a bunch of divides
tomwalters@0: 
tomwalters@0: % limit to keep the poles reasonably stable:
tomwalters@0: interpfactor = (pdamps - mindamp)*damprate; 
tomwalters@0: % interpfactor = min(1.5, interpfactor); % seems to keep it stable
tomwalters@0: 
tomwalters@0: x = xmin + (xmax - xmin).*interpfactor;
tomwalters@0: r = rmin + (rmax - rmin).*interpfactor;
tomwalters@0: 
tomwalters@0: %optional improvement to constellation adds a bit to r:
tomwalters@0: fd = pfreqs.*pdamps;
tomwalters@0: r = r + 0.25*fd.*min(0.05, fd); % quadratic for small values, then linear
tomwalters@0: 
tomwalters@0: zb1 = -2*x;
tomwalters@0: zb2 = r.*r;
tomwalters@0: 
tomwalters@0: bExtraDelay = 1; % boolean to control whether to do allow extra sample of filter delay
tomwalters@0: 
tomwalters@0: if bExtraDelay
tomwalters@0:     % canonic  poles but with input provided where unity DC gain is assured
tomwalters@0:     % (mean value of state is always equal to mean value of input)
tomwalters@0:     newstate = inputs - (state1 - inputs).*zb1 - (state2 - inputs).*zb2;
tomwalters@0:     
tomwalters@0:     % canonic zeros part as before:
tomwalters@0: 
tomwalters@0:     outputs = za0 .* newstate + za1 .* state1 + za2 .* state2 ;
tomwalters@0:     outputs = outputs - 0.0001*outputs.^3; % cubic compression nonlinearity
tomwalters@0: 
tomwalters@0:     state2 = state1;
tomwalters@0:     state1 = newstate;
tomwalters@0: else
tomwalters@0:     input = inputs(1,:); % first stage input is top input; ignore the rest
tomwalters@0:     for stage = 1:nstages
tomwalters@0:         newstate = input - (state1(stage,:) - input).*zb1(stage) - (state2(stage,:) - input).*zb2(stage);
tomwalters@0:         % output of stage immediately becomes input for next stage, if needed:
tomwalters@0:         input = za0(stage) .* newstate + za1(stage) .* state1(stage,:) + za2(stage) .* state2(stage,:) ;
tomwalters@0:         input = input - 0.0001*input.^3; % cubic compression nonlinearity
tomwalters@0:         outputs(stage,:) = input;
tomwalters@0:         state2(stage,:) = state1(stage,:);
tomwalters@0:         state1(stage,:) = newstate;
tomwalters@0:     end
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: if showconstellation
tomwalters@0:     figure(showconstellation);
tomwalters@0:     % keyboard
tomwalters@0:     hold  off
tomwalters@0:     % I have exact x and r, just need y
tomwalters@0:     py = (r.*r - x.*x).^0.5;
tomwalters@0:     % let it error out if we get a complex y.
tomwalters@0:     plot(x,py,'x');
tomwalters@0:     hold on
tomwalters@0:     zr = (za2./za0).^0.5;
tomwalters@0:     zx = -0.5*(za1./za0);
tomwalters@0:     zy = (zr.*zr - zx.*zx).^0.5;
tomwalters@0:     %     zy(zy<0) = 0;
tomwalters@0:     %     zy = zy.^0.5;
tomwalters@0:     plot(zx, zy, 'o');
tomwalters@0:     
tomwalters@0:     drawnow
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if showtransfuns % 1D indexing does only 1 channel for now
tomwalters@0:     fbasis = pi*2.^((-129:0)' /12); %% 120 semitones, 10 octaves, leading up to pi rad/samp
tomwalters@0:     z = exp(i*fbasis);
tomwalters@0:     z2 = z.*z;
tomwalters@0:     nch = size(pfreqs,1);
tomwalters@0:     nf = length(fbasis);
tomwalters@0:     chanxfns = zeros(nf, nch);
tomwalters@0:     cascxfns = ones(nf, nch);
tomwalters@0:     for ch = 1:length(pfreqs)
tomwalters@0:         chanxfns(:,ch) = abs((1+zb1(ch)+zb2(ch)).*(za0(ch)+za1(ch).*z+za2(ch).*z2)./(1+zb1(ch).*z+zb2(ch).*z2));
tomwalters@0:         cascxfns(:,ch) = cascxfns(:,max(1,ch-1)) .* chanxfns(:,ch);
tomwalters@0:     end
tomwalters@0: 
tomwalters@0:     figure(showtransfuns)
tomwalters@0:     hold off
tomwalters@0:     loglog(fbasis, chanxfns);
tomwalters@0:     axis([min(fbasis), max(fbasis), 0.1, 10]);
tomwalters@0: 
tomwalters@0:     figure(showtransfuns+1)
tomwalters@0:     hold off
tomwalters@0:     loglog(fbasis, cascxfns);
tomwalters@0:     axis([min(fbasis), max(fbasis), 0.02, 20000]);
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: return;
tomwalters@0: