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