%% 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
}

You are viewing proxied material from ftp.icm.edu.pl. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.