%% This is the package dspfunctions %% %% (c) Paolo Prandoni %% %% 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 }