path: root/graphics/pstricks/contrib/dsptricks/dspfunctions.sty

%% This is the package dspfunctions
%%
%% (c) Paolo Prandoni <paolo.prandoni _at_ epfl.ch>
%%
%% This program can be redistributed and/or modified under the terms
%% of the LaTeX Project Public License Distributed from CTAN archives
%% in directory macros/latex/base/lppl.txt.
%%
%% DESCRIPTION:
%%   `dspfunctions' is a companion package to dsptricks; it contains a 
%%     set of postscript macros to compute the value of various DSP
%%     common functions
%%
%% v1.1, November 2023
%%

\ProvidesPackage{dspfunctions}[2023/11/08 package for signal processing graphics]

\def\dspToDeg{180 mul }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% \dspRect{a}{b}  rect((x-a)/b)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspRect#1#2{ #1 sub abs #2 div 0.5 gt {0} {1} ifelse }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% \dspTri{a}{b}  triangle((x-a)/b)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspTri#1#2{ #1 sub abs #2 div dup 1 gt {pop 0} {1 exch sub} ifelse }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% \dspExpDec{a}{b}  b^(x-a)u[x-a]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspExpDec#1#2{ #1 sub dup 0 lt {pop 0} {#2 exch exp} ifelse }


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% \dspQuad{a}{b}  quadratic((x-a)/b)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspQuad#1#2{ #1 sub abs #2 div dup 1 gt {pop 0} {dup mul 1 exch sub } ifelse }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Porkpie hat shape (useful for spectral prototypes
% \dspPorkpie{a}{b}  phi((x-a)/b)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspPorkpie#1#2{ #1 sub #2 div dup abs 1 gt {pop 0}%
  {32.4 mul dup cos exch %
            dup 3 mul cos 2 mul exch %
            12 mul cos -0.7 mul  %
            add add 0.31 mul 0.017 add } ifelse }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Raised cosine
% \dspRaisedCos{cutoff}{rolloff}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspRaisedCos#1#2{abs %
  dup dup
    1 #2 sub #1 mul lt %
      {pop pop 1} %
      {1 #2 add #1 mul gt %
        {pop 0} {1 #2 sub #1 mul sub 2 #2 #1 mul mul div 3.14 mul RadtoDeg cos 1 add 0.5 mul} ifelse } %
    ifelse }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Raised cosine, better syntax
% \dspRaisedCos{a}{b}{r} b = cutoff, r = rolloff
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspRaisedCosine#1#2#3{%
  #1 sub abs #2 div
    dup dup
      1 #3 sub lt
        {pop pop 1}
        {1 #3 add gt
          {pop 0}
          {1 #3 sub sub 2 #3 mul div 180 mul cos 1 add 0.5 mul}
          ifelse}
        ifelse }


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% \dspSinc{a}{b}  sinc((x-a)/b)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspSinc#1#2{ #1 sub #2 div dup 0 eq {pop 1} {dup 180 mul sin exch 3.1415 mul div} ifelse }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% \dspSincN{a}{b}  (1/b)sinc((x-a)/b)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspSincN#1#2{\dspSinc{#1}{#2} #2 div }


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% \dspAudio{a}{b}  zero-DC baseband shape (e.g. an audio signal)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspAudio#1#2{% 
  #1 sub #2 div 
  abs 
  dup 1 gt 
    {pop 0}%
    {dup 0.01 le 
      {pop 0}
      {dup 0.2 lt
        {-30 mul}
        {1 exch sub -8 mul}
        ifelse
        2.7 exch exp 1 add 1 exch div 0.5 sub 2 mul
      }
      ifelse}
    ifelse}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Fourier transform of a symmetric 2N+1 tap rect
% \dspSincS{a}{N}  sin((x-a)(2N+1)/2)/sin((x-a)/2)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspSincS#1#2{ #1 sub 90 mul dup #2 2 mul 1 add mul sin exch sin dup 0 eq {pop pop #2 2 mul 1 add} {div} ifelse}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Fourier transform magnitude of a causal N tap rect
% (phase is e^{-j\frac{N-1}{2}\omega})
% \dspSincC{a}{N}  sin((x-a)(N/2))/sin((x-a)/2)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspSincC#1#2{ #1 sub 90 mul dup #2 mul sin exch sin dup 0 eq {pop pop #2} {div} ifelse}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% \dspRand  % Random number uniformly distributed over [-1 1]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspRand{rand 2147483647 div 0.5 sub 2 mul }


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Discrete Fourier Transform of an input sequence; input is  
%  the integer value of the DFT coefficient.
%
% \dspDFTRE{a_0 a_1 ... a_{N-1}}  (real part)
% \dspDFTIM{a_0 a_1 ... a_{N-1}}  (imaginary part)
% \dspDFTMAG{a_0 a_1 ... a_{N-1}} (magnitude)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspDFT#1{%
  cvi [#1] length mod         % DFT is N periodic
  360 mul [#1] length div     % w = (2k\pi)/N
  0                           % index n
  0                           % accumulator Re
  0                           % accumulator Im
  [#1]                        % data points
  {                           %  STACK:
                              % w n re im a_n
    dup                       % w n re im a_n a_n
    5 index                   % w n re im a_n a_n w
    5 index                   % w n re im a_n a_n w n
    mul dup                   % w n re im a_n a_n nw nw
    sin exch cos              % w n re im a_n a_n sin(nw) cos(nw)
    4 1 roll mul              % w n re im cos(nw) a_n (a_n)sin(nw)
    3 1 roll mul              % w n re im (a_n)sin(nw) (a_n)cos(nw)
    4 1 roll add              % w n (a_n)cos(nw) re im'
    3 1 roll add exch         % w n re' im'
    3 2 roll 1 add            % w re' im' n'
    3 1 roll                  % w n re im
  } forall
  4 2 roll pop pop            % re im
}
\def\dspDFTRE#1{\dspDFT{#1} pop }
\def\dspDFTIM#1{\dspDFT{#1} exch pop }
\def\dspDFTMAG#1{\dspDFT{#1} dup mul exch dup mul add sqrt }


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Frequency response of a (2N+1)-tap Type-I FIR computed at a given 
%  normalized frequency value. Frequency response is real for Type-I
% The filter is considered zero-centered, so a_0 is the center tap
%
% \dspFIRI{a_0 a_1 ... a_{N-1}}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspFIRI#1{%
  \dspToDeg                   % input to degrees  
  0                           % index n
  0                           % accumulator A
  [#1]                        % coefficients a_n
  {
    3 index                   % x
    3 index                   % n [*** using index INCREASES stack size... so it's 3 3 rather than 3 2]
    mul cos mul               % a_n cos nx
    add                       % accumulate
    exch 1 add exch           % i++
  } forall
  3 1 roll pop pop
  2 mul                       % final value is 2A - a_0
  [#1] 0 get sub
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Magnitude response of a generic digital filter defined by the 
%  constant-coefficient difference equation:
%   y[n] = b_0 x[n] + b_1 x[n-1] + ... + b_{N-1} x[n-N+1] 
%          - a_1 y[n-1] - ... - a_{M-1} y[n-M+1]
%
% The response is computed at the given normalized frequency value
%
% \dspTFM{b_0 b_1 b_2 ... b_{M-1}}{a_0 a_1 ... a_{N-1}}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspTFM#1#2{%
  \dspToDeg                   % input to degrees
  dup                         % save a copy for denominator
  0                           % index n
  0                           % accumulator Re
  0                           % accumulator Im
  [#1]                        % coefficients b_n
  {                           %  STACK (neglecting saved input at bottom):
                              % x n re im b_n
    dup                       % x n re im b_n b_n
    5 index                   % x n re im b_n b_n x
    5 index                   % x n re im b_n b_n x n
    mul dup                   % x n re im b_n b_n nx nx
    sin exch cos              % x n re im b_n b_n sin(nx) cos(nx)
    4 1 roll mul              % x n re im cos(nx) b_n (b_n)sin(nx)
    3 1 roll mul              % x n re im (b_n)sin(nx) (b_n)cos(nx)
    4 1 roll add              % x n (b_n)cos(nx) re im'
    3 1 roll add exch         % x n re' im'
    3 2 roll 1 add            % x re' im' n'
    3 1 roll                  % x n re im
  } forall
  4 2 roll pop pop            % re im
  dup mul exch dup mul add    % (re^2 + im^2)
  sqrt                        % mag of the numerator of transfer function
  exch                        % bring up saved input copy
  0                           % same loop for the a_n coefficients
  0                           
  0                           
  [1 #2]                        
  {                           
    dup                       
    5 index                   
    5 index                   
    mul dup                  
    sin exch cos              
    4 1 roll mul              
    3 1 roll mul              
    4 1 roll add              
    3 1 roll add exch         
    3 2 roll 1 add            
    3 1 roll                  
  } forall
  4 2 roll pop pop
  dup mul exch dup mul add 
  sqrt                            
  div %0 eq {pop pop 0} {div} ifelse
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Filter the data with a system implementing the transfer function
%   y[n] = b_0 x[n] + b_1 x[n-1] + ... + b_{N-1} x[n-N+1] 
%          - a_1 y[n-1] - ... - a_{M-1} y[n-M+1]
%
% Use the setFilter macro in the setup part of the drawing command
%  and then \dspFilter{b_0 b_1 b_2 ... b_{M-1}}{a_1 ... a_{N-1}}
%
% \dspSignalOpt{\dspSetFilter{b_0 b_1 b_2 ... b_{M-1}}{a_1 ... a_{N-1}}}{x ... \dspFilter}
%
% NB: skip a_0 (assumed = 1) from the feedback part!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\dspSetFilter#1#2{%
  /b [#1] def 
  /xlen b length def 
  /xbuf [xlen {0} repeat] def
  /a [#2] def 
  /ylen a length def 
  /ybuf [ylen {0} repeat] def
}
\def\dspFilter{
  xbuf aload pop pop xbuf astore /xbuf exch def 
  0
  0 1 xlen 1 sub {
    dup
    xbuf exch get
    exch
    b exch get
    mul add
  } for
  0 1 ylen 1 sub {
    dup
    ybuf exch get
    exch
    a exch get
    mul sub
  } for
  dup
  ybuf aload pop pop ybuf astore /ybuf exch def 
}