blob: ae63eccaa0dc39640364b1ab2c23f4005ea5f64c (
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
|
%% 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.0, July 2014
%%
\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 }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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] = a_0 x[n] + a_1 x[n-1] + ... + a_{N-1} x[n-N+1]
% - b_1 y[n-1] - ... - b_{M-1} y[n-M+1]
%
% The response is computed at the given normalized frequency value
%
% \dspTFM{a_0 a_1 ... a_{N-1}}{b_1 b_2 ... b_{M-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 a_n
{ % STACK (neglecting saved input at bottom):
% x n re im a_n
dup % x n re im a_n a_n
5 index % x n re im a_n a_n x
5 index % x n re im a_n a_n x n
mul dup % x n re im a_n a_n nx nx
sin exch cos % x n re im a_n a_n sin(nx) cos(nx)
4 1 roll mul % x n re im cos(nx) a_n (a_n)sin(nx)
3 1 roll mul % x n re im (a_n)sin(nx) (a_n)cos(nx)
4 1 roll add % x n (a_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 b_n coefficients
0
0
[#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
}
|
