function [y,Zf,Zmax]=resfiltq(wm,rm,zm,x,Zi,Q,runmod)
%RESFILTQ  Resonator-based filtering with quantization
%
% See also: tf2res, res2tf
%
%Input parameters
%   wm - weighting factors, as provided by tf2res
%   rm - weighting factors, as provided by tf2res
%   zm - resonator poles, as provided by tf2res
%   x - input vector
%   Zi - Zi output of the previous run (as in filter)
%   Q - quantizer for the simulation: see roundoff tb (double for empty Q)
%   runmod - run modifier structure. Fields:
%     implementation: 'real', 'direct II' (same: real resonators)
%         'coupled': resonators calculated as complex pairs
%         'complex': resonators calculated as complex resonators
%     Qmul - quantizer after multiplications
%     Qadd - quantizer after additions
%     Qcoeff - quantizer of the constants

%	Copyright (c) I. Kollar, 2008, all rights reserved
%	Last modified: 20-Apr-2010

%Literature:
%Gábor Péceli, "Resonator-based Digital Filters," IEEE Trans. on
%Circuits and Systems, Vol. 36, No. 1, Jan. 1989, pp. 156-159.

global ovf_global
ovf_global_save=ovf_global; ovf_global=[];
global resfiltq_sign_of_life

error(nargchk(4,7,nargin))
N=length(wm);
if isempty(rm), rm=1/N*ones(N,1); end
if length(wm)==length(rm)+1, b0=wm(end); wm(end)=[]; N=length(wm); else b0=[]; end
if isempty(b0), b0=0; end

if length(rm)~=length(wm), error('lengths of rm and wm differ'), end
if isempty(zm), zm=exp(j*[0:N-1]'/N*2*pi); end
if length(zm)~=length(wm), error('lengths of zm and wm differ'), end

if nargin<5, Zi=[]; end, if isempty(Zi), Zi=zeros(size(wm)); end
if nargin<6, Q=[]; end
if nargin<7
  if (nargin>=6)&isstruct(Q)&isfield(Q,'implementation'), runmod=Q; Q=[];
  else runmod=[];
  end
end

if ~isfield(runmod,'implementation'), runmod.implementation='real'; end
if ~isfield(runmod,'Qmul'), runmod.Qmul=[]; end
Qmul=runmod.Qmul; if isempty(Qmul), Qmul=Q; end
if ~isfield(runmod,'Qadd'), runmod.Qadd=[]; end
Qadd=runmod.Qadd; if isempty(Qadd), Qadd=Q; end
if ~isfield(runmod,'Qcoeff'), runmod.Qcoeff=[]; end
Qcoeff=runmod.Qcoeff; if isempty(Qcoeff), Qcoeff=Q; end
if ~qstructcmp(Qcoeff,'double')
  %wm=roundq(wm,Qcoeff,'complex');
  ind=find(~ismember([0:N-1],N/4*[0:4]));
  zm(ind)=roundq(zm(ind),Qcoeff,'complex'); %do not destroy zm=1
end
if any(ovf_global), warning('Overflow at initializing'), end
if nargout>=3, Zmax=[inf,-inf]; end
y=zeros(size(x));
if rem(N,2)==0
  zm(N/2+1)=real(zm(N/2+1)); wm(N/2+1)=real(wm(N/2+1));
  rm(N/2+1)=real(rm(N/2+1));
end
if strcmpi(runmod.implementation,'real')|strncmpi(runmod.implementation,'direct II',3)|isempty(runmod.implementation)
  %rm=real(gm./zm);
  phirm=angle(rm); rm=abs(rm);
  if ~qstructcmp(Qcoeff,'double')
    rm=roundq(rm,Qcoeff);
    %phirm=roundq(phirm,Qcoeff);
  end
  for n=1:length(y)
    %Feedback and output calculation
    fb=0; %y(n)=0;
    for m=0:N/2
      if (m==0)|(m==N/2) %1st order
        fb=roundq(fb+Zi(m+1),Qadd);
        y(n)=roundq(y(n)+roundq(wm(m+1)*Zi(m+1),Qmul),Qadd);
      else %2nd order
        m2=N-m; cosphim=real(zm(m+1)); sinphim=imag(zm(m+1));
        %2*Re x
        if all(phirm==0) %&0
          fb=roundq(fb+roundq(2*roundq(roundq(cosphim*Zi(m+1),Qmul)...
            -Zi(m2+1),Qadd),Qmul),Qadd);
          %2*Re x*w
          y(n)=roundq(y(n)+roundq(2*...
            roundq(...
            roundq(roundq(roundq(cosphim*Zi(m+1),Qmul)...
            -Zi(m2+1),Qadd)*real(wm(m+1)),Qmul)...
            -roundq(roundq(sinphim*Zi(m2+1),Qmul)*imag(wm(m+1)),Qmul),...
            Qadd),Qmul),Qadd);
        else %complex rm
          Rex=roundq(roundq(cosphim*Zi(m+1),Qmul)-Zi(m2+1),Qadd);
          Imx=roundq(sinphim*Zi(m+1),Qmul);
          fb=roundq(fb+roundq(2*real(cmulq(exp(j*phirm(m+1)),Rex+j*Imx,Qadd)),Qmul),Qadd);
          %2*Re x*w
          y(n)=roundq(y(n)+roundq(2*...
            real(cmulq(cmulq(exp(j*phirm(m+1)),Rex+j*Imx,Qadd),wm(m+1),Qadd)),...
            Qmul),Qadd);
        end
      end %rm real or complex
    end
    %fb=0; %for testing only
    fin=roundq(x(n)-fb,Qadd);
    y(n)=roundq(y(n)+roundq(b0*x(n),Qmul),Qadd);
    if ~isequal(Qadd,Q)&~qstructcmp(Qadd,Q), y(n)=roundq(y(n),Q); end
    %new Zi's
    for m=0:N/2
      if (m==0)|(m==N/2) %1st order
        Zi(m+1)=roundq(roundq(rm(m+1)*zm(m+1)*fin,Qmul)+...
          roundq(zm(m+1)*Zi(m+1),Qmul),Qadd);
        %For zm=1, this is OK for fixed-point
      else %2nd order
        m2=N-m; cosphim=real(zm(m+1)); sinphim=imag(zm(m+1));
        Zim1=roundq(roundq(roundq(2*roundq(Zi(m+1)*cosphim,Qmul),Qmul)-Zi(m2+1))...
          +roundq(rm(m+1)*fin,Qmul),Qadd);
        Zi(m2+1)=Zi(m+1);
        if ~isequal(Qadd,Q)&~qstructcmp(Qadd,Q), Zim1=roundq(Zim1,Q); end
        Zi(m+1)=Zim1;
      end
    end
    if ~isequal(Qadd,Q)&~qstructcmp(Qadd,Q), Zi=roundq(Zi,Q); end
    if nargout>=3
      Zmax(2)=max(Zmax(2),max(Zi));
      Zmax(1)=min(Zmax(1),min(Zi));
    end
    sig='.'; if any(ovf_global), sig='!'; end
    if any(resfiltq_sign_of_life)&(rem(n,20)==0), fprintf(sig), end
  end %for n
elseif strcmpi(runmod.implementation,'coupled')
  %Gábor Péceli, "Structurally Passive Resonator-based Digital Filters,"
  %ISCAS'88, 1988 International Symposium on Circuits and Systems,
  %pp. 1421-1424, Fig. 2b.
  %error('Coupled not yet ready')
  phirm=angle(rm); rm=abs(rm);
  if ~qstructcmp(Qcoeff,'double')
    rm=roundq(rm,Qcoeff);
    %phirm=roundq(phirm,Qcoeff);
  end
  for n=1:length(y)
    %Feedback and output calculation
    fb=0; %y(n)=0;
    for m=0:N/2
      if (m==0)|(m==N/2) %1st order
        fb=roundq(fb+Zi(m+1),Qadd);
        y(n)=roundq(y(n)+roundq(wm(m+1)*Zi(m+1),Qmul),Qadd);
      else %2nd order
        m2=N-m; cosphim=real(zm(m+1)); sinphim=imag(zm(m+1));
        %2*Re x
        if phirm(m+1)==0 %real rm
          fb=roundq(fb+roundq(2*roundq(roundq(cosphim*Zi(m+1),Qmul)...
            +roundq(sinphim*Zi(m2+1),Qmul),Qadd),Qmul),Qadd);
          %2*Re x*w
          y(n)=roundq(y(n)+roundq(2*roundq(...
            roundq(roundq(roundq(cosphim*Zi(m+1),Qmul)...
            +roundq(sinphim*Zi(m2+1),Qmul),Qadd)*real(wm(m+1)),Qmul)...
            +roundq(roundq(roundq(-sinphim*Zi(m+1),Qmul)...
            +roundq(cosphim*Zi(m2+1),Qmul),Qadd)*imag(wm(m+1)),Qmul),...
            Qadd),Qmul),Qadd);
        else %complex rm
          Rex=roundq(roundq(cosphim*Zi(m+1),Qmul)...
            +roundq(sinphim*Zi(m2+1),Qmul),Qadd);
          Imx=-roundq(roundq(-sinphim*Zi(m+1),Qmul)...
            +roundq(cosphim*Zi(m2+1),Qmul),Qadd);
          fb=roundq(fb+roundq(2*real(cmulq(exp(j*phirm(m+1)),Rex+j*Imx,Qadd)),...
            Qmul),Qadd);
          %2*Re x*w
          y(n)=roundq(y(n)+roundq(2*real(cmulq(cmulq(exp(j*phirm(m+1)),Rex+j*Imx,Qadd),...
            wm(m+1),Qadd)),Qmul),Qadd);
        end
      end %rm real or complex
    end
    %fb=0; %for testing only
    fin=roundq(x(n)-fb,Qadd);
    y(n)=roundq(y(n)+roundq(b0*x(n),Qmul),Qadd);
    if ~isequal(Qadd,Q)&~qstructcmp(Qadd,Q), y(n)=roundq(y(n),Q); end
    %new Zi's
    for m=0:N/2
      if (m==0)|(m==N/2) %1st order
        Zi(m+1)=roundq(roundq(rm(m+1)*zm(m+1)*fin,Qmul)+...
          roundq(zm(m+1)*Zi(m+1),Qmul),Qadd);
        %For zm=1, this is OK for fixed-point
      else %2nd order
        m2=N-m; cosphim=real(zm(m+1)); sinphim=imag(zm(m+1));
        Zim1=roundq(roundq(roundq(Zi(m+1)*cosphim,Qmul)...
          +roundq(Zi(m2+1)*sinphim,Qmul),Qadd)...
          +roundq(rm(m+1)*fin,Qmul),Qadd);
        Zim2=roundq(roundq(roundq(Zi(m2+1)*cosphim,Qmul)...
          +roundq(Zi(m+1)*(-sinphim),Qmul),Qadd),Qadd);
        if ~isequal(Qadd,Q)&~qstructcmp(Qadd,Q), Zim1=roundq(Zim1,Q); end
        Zi(m+1)=Zim1;
        if ~isequal(Qadd,Q)&~qstructcmp(Qadd,Q), Zim2=roundq(Zim2,Q); end
        Zi(m2+1)=Zim2;
      end
    end %for m
    if ~isequal(Qadd,Q)&~qstructcmp(Qadd,Q), Zi=roundq(Zi,Q); end
    if nargout>=3
      Zmax(2)=max(Zmax(2),max(Zi));
      Zmax(1)=min(Zmax(1),min(Zi));
    end
    sig='.'; if any(ovf_global), sig='!'; end
    if any(resfiltq_sign_of_life)&(rem(n,20)==0), fprintf(sig), end
  end %for n
elseif strcmpi(runmod.implementation,'complex')
  if ~qstructcmp(Qcoeff,'double')
    rm=roundq(rm,Qcoeff);
  end
  for n=1:length(y)
    %Feedback and output calculation
    fb=0; %y(n)=0;
    for m=0:N/2
      if (m==0)|(m==N/2) %1st order
        fb=roundq(fb+Zi(m+1),Qadd);
        if isfield(Qmul,'Type')&strncmpi(Qmul.Type,'fixed',2)
          mx=Qmul.Maximum;
          if abs(real(wm(m+1)))<=mx
            y(n)=roundq(y(n)+roundq(wm(m+1)*Zi(m+1),Qmul),Qadd);
          elseif abs(real(wm(m+1)))<=2*mx %try to execute with w's over 1
            term=roundq(wm(m+1)/2*Zi(m+1),Qmul);
            y(n)=roundq(roundq(y(n)+term,Qadd)+term,Qadd);
          else error('wm overflow')
          end
        else
          y(n)=roundq(y(n)+roundq(wm(m+1)*Zi(m+1),Qmul),Qadd);
        end
      else %complex conjugate pair
        fb=roundq(fb+roundq(2*real(Zi(m+1)),Qmul,'complex'),Qadd,'complex');
        y(n)=roundq(y(n)+roundq(2*real(roundq(wm(m+1)*Zi(m+1),...
          Qmul,'complex')),Qmul,'complex'),Qadd);
      end
    end
    %fb=0; %for testing only
    fin=roundq(x(n)-fb,Qadd);
    y(n)=roundq(y(n)+roundq(b0*x(n),Qmul),Qadd);
    %new Zi's
    for m=0:N/2
      if ismember(m,[0,N/2])
        Zi(m+1)=roundq(roundq(rm(m+1)*zm(m+1)*fin,Qmul,'complex')+...
          roundq(zm(m+1)*Zi(m+1),Qmul,'complex'),Qadd,'complex');
        %For zm=1, this is OK for fixed-point
      else %complex conjugate pair
        Zi(m+1)=roundq(roundq(roundq(rm(m+1)*zm(m+1),Qmul,'complex')*fin,Qmul,'complex')...
          +roundq(zm(m+1)*Zi(m+1),Qmul,'complex'),Qadd,'complex');
      end
    end
    Zi(N:-1:ceil(N/2)+2)=conj(Zi(2:floor(N/2)));
    if ~isequal(Qadd,Q)&~qstructcmp(Qadd,Q)
      Zi=roundq(Zi,Q,'complex');
      y(n)=roundq(y(n),Q);
    end
    if nargout>=3
      Zmax(2)=max([Zmax(2),max(real(Zi)),max(imag(Zi))]);
      Zmax(1)=min([Zmax(1),min(real(Zi)),min(imag(Zi))]);
    end
    sig='.'; if any(ovf_global), sig='!'; end
    if any(resfiltq_sign_of_life)&(rem(n,20)==0), fprintf(sig), end
  end %for n
else
  error('runmod.implementation is not allowed')
end
if any(resfiltq_sign_of_life), fprintf('\n'), end
if nargout>=2, Zf=Zi; end

if any(ovf_global), warning('Overflow happened in resfiltq'), end
if any(ovf_global_save), ovf_global=ovf_global_save; end

Error using ==> resfiltq at 32
Not enough input arguments.

Published with MATLAB® 7.10