summaryrefslogtreecommitdiff
path: root/graphics/pstricks/contrib/dsptricks/dspfunctions.sty
blob: c42588a55d15a52c3b92208dc9a7d8ef534ddae3 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
%% 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}[2014/07/24 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 
}