1 |
root |
1.8 |
$VERSION = '1.1'; |
2 |
root |
1.1 |
|
3 |
|
|
pp_setversion $VERSION; |
4 |
|
|
pp_beginwrap (); # force error with older PPs |
5 |
|
|
|
6 |
|
|
pp_addpm {At => Top}, <<'EOD'; |
7 |
|
|
=head1 NAME |
8 |
|
|
|
9 |
|
|
PDL::Audio - Some PDL functions intended for audio processing. |
10 |
|
|
|
11 |
|
|
=head1 SYNOPSIS |
12 |
|
|
|
13 |
|
|
use PDL; |
14 |
|
|
use PDL::Audio; |
15 |
|
|
|
16 |
|
|
=head1 DESCRIPTION |
17 |
|
|
|
18 |
root |
1.6 |
Oh well ;) Not much "introductory documentation" has been written yet :( |
19 |
root |
1.1 |
|
20 |
|
|
=head2 NOTATION |
21 |
|
|
|
22 |
|
|
Brackets around parameters indicate that the respective parameter is |
23 |
|
|
optional and will be replaced with some default value when absent (or |
24 |
|
|
C<undef>, which might be different in other packages). |
25 |
|
|
|
26 |
|
|
The sampling frequency and duration are by default (see individual |
27 |
|
|
descriptions) given in cycles/sample (or samples in case of a |
28 |
|
|
duration). That means if you want to specify a duration of two seconds, |
29 |
|
|
you have to multiply by the sampling frequency in HZ, and if you want |
30 |
|
|
to specify a frequency of 440 Hz, you have to divide by the sampling |
31 |
|
|
frequency: |
32 |
|
|
|
33 |
|
|
# Syntax: gen_oscil duration*, frequency/ |
34 |
|
|
$signal = gen_oscil 2*HZ, 440/HZ; |
35 |
|
|
# with a sampling frequency of 44100 Hertz: |
36 |
|
|
$signal = gen_oscil 2*44100, 440/44100; |
37 |
|
|
|
38 |
|
|
To help you, the required unit is given as a type suffix in the parameter |
39 |
|
|
name. A "/" means that you have to divide by the sampling frequency (to |
40 |
|
|
convert from Hertz) and a suffix of "*" indicates that a multiplication is |
41 |
|
|
required. |
42 |
|
|
|
43 |
|
|
Most parameters named "size", "duration" (or marked with "*") can |
44 |
|
|
be replaced by a piddle, which is then used to give length and from |
45 |
|
|
(mono/stereo). |
46 |
|
|
|
47 |
|
|
=head2 HEADER ATTRIBUTES |
48 |
|
|
|
49 |
|
|
The following header attributes are stored and evaluated by most |
50 |
|
|
functions. PDL::Audio provides mutator methods for all them (e.g. |
51 |
|
|
|
52 |
|
|
print "samplerate is ", $pdl->rate; |
53 |
|
|
$pdl->comment("set the comment to this string"); |
54 |
|
|
|
55 |
|
|
=over 4 |
56 |
|
|
|
57 |
|
|
=item rate |
58 |
|
|
|
59 |
|
|
The sampling rate in hz. |
60 |
|
|
|
61 |
|
|
=item filetype |
62 |
|
|
|
63 |
|
|
The filetype (wav, au etc..). Must be one of: |
64 |
|
|
|
65 |
|
|
FILE_NEXT FILE_AIFC FILE_RIFF FILE_BICSF FILE_NIST FILE_INRS FILE_ESPS |
66 |
|
|
FILE_SVX FILE_VOC FILE_SNDT FILE_RAW FILE_SMP FILE_SD2 FILE_AVR |
67 |
|
|
FILE_IRCAM FILE_SD1 FILE_SPPACK FILE_MUS10 FILE_HCOM FILE_PSION |
68 |
|
|
FILE_MAUD FILE_IEEE FILE_DESKMATE FILE_DESKMATE_2500 FILE_MATLAB |
69 |
|
|
FILE_ADC FILE_SOUNDEDIT FILE_SOUNDEDIT_16 FILE_DVSM FILE_MIDI |
70 |
|
|
FILE_ESIGNAL FILE_SOUNDFONT FILE_GRAVIS FILE_COMDISCO FILE_GOLDWAVE |
71 |
|
|
FILE_SRFS FILE_MIDI_SAMPLE_DUMP FILE_DIAMONDWARE FILE_REALAUDIO |
72 |
|
|
FILE_ADF FILE_SBSTUDIOII FILE_DELUSION FILE_FARANDOLE FILE_SAMPLE_DUMP |
73 |
|
|
FILE_ULTRATRACKER FILE_YAMAHA_SY85 FILE_YAMAHA_TX16 FILE_DIGIPLAYER |
74 |
|
|
FILE_COVOX FILE_SPL FILE_AVI FILE_OMF FILE_QUICKTIME FILE_ASF |
75 |
|
|
FILE_YAMAHA_SY99 FILE_KURZWEIL_2000 FILE_AIFF FILE_AU |
76 |
|
|
|
77 |
|
|
=item path |
78 |
|
|
|
79 |
|
|
The filename (or file specification) used to load or save a file. |
80 |
|
|
|
81 |
|
|
=item format |
82 |
|
|
|
83 |
|
|
Specifies the type the underlying file format uses. The samples will always be in short or long signed format. |
84 |
|
|
|
85 |
|
|
Must be one of |
86 |
|
|
|
87 |
|
|
FORMAT_NO_SND FORMAT_16_LINEAR FORMAT_8_MULAW FORMAT_8_LINEAR |
88 |
|
|
FORMAT_32_FLOAT FORMAT_32_LINEAR FORMAT_8_ALAW FORMAT_8_UNSIGNED |
89 |
|
|
FORMAT_24_LINEAR FORMAT_64_DOUBLE FORMAT_16_LINEAR_LITTLE_ENDIAN |
90 |
|
|
FORMAT_32_LINEAR_LITTLE_ENDIAN FORMAT_32_FLOAT_LITTLE_ENDIAN |
91 |
|
|
FORMAT_64_DOUBLE_LITTLE_ENDIAN FORMAT_16_UNSIGNED |
92 |
|
|
FORMAT_16_UNSIGNED_LITTLE_ENDIAN FORMAT_24_LINEAR_LITTLE_ENDIAN |
93 |
|
|
FORMAT_32_VAX_FLOAT FORMAT_12_LINEAR FORMAT_12_LINEAR_LITTLE_ENDIAN |
94 |
|
|
FORMAT_12_UNSIGNED FORMAT_12_UNSIGNED_LITTLE_ENDIAN COMPATIBLE_FORMAT |
95 |
|
|
|
96 |
|
|
PDL::Audio conviniently defines the following aliases for the following |
97 |
|
|
constants, that are already correct for the host byteorder: |
98 |
|
|
|
99 |
|
|
FORMAT_ULAW_BYTE FORMAT_ALAW_BYTE FORMAT_LINEAR_BYTE |
100 |
|
|
FORMAT_LINEAR_SHORT FORMAT_LINEAR_USHORT FORMAT_LINEAR_LONG |
101 |
|
|
FORMAT_LINEAR_FLOAT FORMAT_LINEAR_DOUBLE |
102 |
|
|
|
103 |
|
|
=item comment |
104 |
|
|
|
105 |
|
|
The file comment (if any). |
106 |
|
|
|
107 |
|
|
=item device |
108 |
|
|
|
109 |
|
|
The device to output audio. One of: |
110 |
|
|
|
111 |
|
|
DEV_DEFAULT DEV_READ_WRITE DEV_ADAT_IN DEV_AES_IN DEV_LINE_OUT |
112 |
|
|
DEV_LINE_IN DEV_MICROPHONE DEV_SPEAKERS DEV_DIGITAL_IN DEV_DIGITAL_OUT |
113 |
|
|
DEV_DAC_OUT DEV_ADAT_OUT DEV_AES_OUT DEV_DAC_FILTER DEV_MIXER |
114 |
|
|
DEV_LINE1 DEV_LINE2 DEV_LINE3 DEV_AUX_INPUT DEV_CD_IN DEV_AUX_OUTPUT |
115 |
|
|
DEV_SPDIF_IN DEV_SPDIF_OUT |
116 |
|
|
|
117 |
|
|
=back |
118 |
|
|
|
119 |
|
|
=head2 EXPORTED CONSTANTS |
120 |
|
|
|
121 |
|
|
In addition to the exported constants described above (and later in |
122 |
|
|
the function descriptions), this module also exports the mathematical |
123 |
|
|
constants M_PI and M_2PI, so watch out for clashes! |
124 |
|
|
|
125 |
|
|
=cut |
126 |
|
|
EOD |
127 |
|
|
|
128 |
|
|
@cconsts = qw( |
129 |
|
|
SNDLIB_DEFAULT_DEVICE SNDLIB_READ_WRITE_DEVICE |
130 |
|
|
SNDLIB_ADAT_IN_DEVICE SNDLIB_AES_IN_DEVICE |
131 |
|
|
SNDLIB_LINE_OUT_DEVICE SNDLIB_LINE_IN_DEVICE |
132 |
|
|
SNDLIB_MICROPHONE_DEVICE SNDLIB_SPEAKERS_DEVICE |
133 |
|
|
SNDLIB_DIGITAL_IN_DEVICE SNDLIB_DIGITAL_OUT_DEVICE |
134 |
|
|
SNDLIB_DAC_OUT_DEVICE SNDLIB_ADAT_OUT_DEVICE |
135 |
|
|
SNDLIB_AES_OUT_DEVICE SNDLIB_DAC_FILTER_DEVICE |
136 |
|
|
SNDLIB_MIXER_DEVICE SNDLIB_LINE1_DEVICE |
137 |
|
|
SNDLIB_LINE2_DEVICE SNDLIB_LINE3_DEVICE |
138 |
|
|
SNDLIB_AUX_INPUT_DEVICE SNDLIB_CD_IN_DEVICE |
139 |
|
|
SNDLIB_AUX_OUTPUT_DEVICE SNDLIB_SPDIF_IN_DEVICE |
140 |
|
|
SNDLIB_SPDIF_OUT_DEVICE |
141 |
|
|
|
142 |
|
|
SNDLIB_NO_SND SNDLIB_16_LINEAR SNDLIB_8_MULAW SNDLIB_8_LINEAR |
143 |
|
|
SNDLIB_32_FLOAT SNDLIB_32_LINEAR SNDLIB_8_ALAW SNDLIB_8_UNSIGNED |
144 |
|
|
SNDLIB_24_LINEAR SNDLIB_64_DOUBLE SNDLIB_16_LINEAR_LITTLE_ENDIAN |
145 |
|
|
SNDLIB_32_LINEAR_LITTLE_ENDIAN SNDLIB_32_FLOAT_LITTLE_ENDIAN |
146 |
|
|
SNDLIB_64_DOUBLE_LITTLE_ENDIAN SNDLIB_16_UNSIGNED |
147 |
|
|
SNDLIB_16_UNSIGNED_LITTLE_ENDIAN SNDLIB_24_LINEAR_LITTLE_ENDIAN |
148 |
|
|
SNDLIB_32_VAX_FLOAT SNDLIB_12_LINEAR SNDLIB_12_LINEAR_LITTLE_ENDIAN |
149 |
|
|
SNDLIB_12_UNSIGNED SNDLIB_12_UNSIGNED_LITTLE_ENDIAN |
150 |
|
|
SNDLIB_COMPATIBLE_FORMAT |
151 |
|
|
|
152 |
|
|
NeXT_sound_file AIFC_sound_file RIFF_sound_file |
153 |
|
|
BICSF_sound_file NIST_sound_file INRS_sound_file |
154 |
|
|
ESPS_sound_file SVX_sound_file VOC_sound_file |
155 |
|
|
SNDT_sound_file raw_sound_file SMP_sound_file |
156 |
|
|
SD2_sound_file AVR_sound_file IRCAM_sound_file |
157 |
|
|
SD1_sound_file SPPACK_sound_file MUS10_sound_file |
158 |
|
|
HCOM_sound_file PSION_sound_file MAUD_sound_file |
159 |
|
|
IEEE_sound_file DeskMate_sound_file DeskMate_2500_sound_file |
160 |
|
|
Matlab_sound_file ADC_sound_file SoundEdit_sound_file |
161 |
|
|
SoundEdit_16_sound_file DVSM_sound_file MIDI_file |
162 |
|
|
Esignal_file soundfont_sound_file gravis_sound_file |
163 |
|
|
comdisco_sound_file goldwave_sound_file srfs_sound_file |
164 |
|
|
MIDI_sample_dump DiamondWare_sound_file RealAudio_sound_file |
165 |
|
|
ADF_sound_file SBStudioII_sound_file Delusion_sound_file |
166 |
|
|
Farandole_sound_file Sample_dump_sound_file Ultratracker_sound_file |
167 |
|
|
Yamaha_SY85_sound_file Yamaha_TX16_sound_file digiplayer_sound_file |
168 |
|
|
Covox_sound_file SPL_sound_file AVI_sound_file |
169 |
|
|
OMF_sound_file Quicktime_sound_file asf_sound_file |
170 |
|
|
Yamaha_SY99_sound_file Kurzweil_2000_sound_file AIFF_sound_file |
171 |
|
|
|
172 |
|
|
BANDPASS DIFFERENTIATOR HILBERT |
173 |
|
|
); |
174 |
|
|
|
175 |
|
|
for (@cconsts) { |
176 |
|
|
local $_ = $_; |
177 |
|
|
s/^SNDLIB_(.*)_DEVICE/DEV_$1/; |
178 |
|
|
s/^SNDLIB_/FORMAT_/; |
179 |
|
|
s/MIDI_file/FILE_MIDI/; |
180 |
|
|
s/MIDI_sample_dump/FILE_MIDI_sample_dump/; |
181 |
|
|
s/Esignal_file/FILE_Esignal/; |
182 |
|
|
s/(.*)_sound_file/FILE_$1/; |
183 |
|
|
push @pconsts, uc $_; |
184 |
|
|
} |
185 |
|
|
|
186 |
|
|
@METHODS = qw( |
187 |
|
|
waudio cut_leading_silence cut_trailing_silence cut_silence filter_zpolar |
188 |
|
|
partials2polynomial scale2short filter_fir filter_iir rfft irfft spectrum |
189 |
|
|
filter_comb filter_notch filter_allpass filter_src filter_granulate |
190 |
|
|
filter_contrast_enhance describe_audio filter_center ring_modulate |
191 |
|
|
amplitude_modulate filter_ppolar design_remez_fir concat |
192 |
|
|
gen_oscil gen_sawtooth gen_square gen_triangle gen_asymmetric_fm gen_rand_1f |
193 |
|
|
gen_env audiomix gen_sum_of_cosines gen_sine_summation gen_fft_window |
194 |
|
|
gen_pulse_train gen_rand gen_from_table gen_from_partials partials2waveshape |
195 |
|
|
|
196 |
|
|
); |
197 |
|
|
@EXPORT = (qw( |
198 |
|
|
sound_format_name sound_type_name raudio |
199 |
|
|
FORMAT_ULAW_BYTE FORMAT_ALAW_BYTE FORMAT_LINEAR_BYTE |
200 |
|
|
FORMAT_LINEAR_SHORT FORMAT_LINEAR_USHORT FORMAT_LINEAR_LONG |
201 |
|
|
FORMAT_LINEAR_FLOAT FORMAT_LINEAR_DOUBLE |
202 |
|
|
M_PI M_2PI |
203 |
|
|
), @METHODS, @pconsts); |
204 |
|
|
|
205 |
|
|
for (@EXPORT) { |
206 |
|
|
pp_add_exported "", $_; |
207 |
|
|
} |
208 |
|
|
|
209 |
|
|
pp_addpm "\@METHODS = qw(@METHODS);\n"; |
210 |
|
|
|
211 |
|
|
$addboot = "{\nHV *stash = gv_stashpvn(\"PDL::Audio\", 10, TRUE);"; |
212 |
|
|
|
213 |
|
|
for (0..$#cconsts) { |
214 |
|
|
$addboot .= " newCONSTSUB(stash, \"$pconsts[$_]\", newSViv($cconsts[$_]));\n"; |
215 |
|
|
} |
216 |
|
|
$addboot .= " newCONSTSUB(stash, \"M_PI\" , newSVnv(M_PI ));\n"; |
217 |
|
|
$addboot .= " newCONSTSUB(stash, \"M_2PI\", newSVnv(M_2PI));\n"; |
218 |
|
|
$addboot .= '} |
219 |
|
|
mus_set_error_handler (mus_barfer); |
220 |
|
|
'; |
221 |
|
|
|
222 |
|
|
pp_add_boot $addboot; |
223 |
|
|
|
224 |
|
|
pp_addpm <<'EOD'; |
225 |
|
|
|
226 |
|
|
use PDL::LiteF; |
227 |
|
|
use PDL::Complex; |
228 |
|
|
use Carp; |
229 |
|
|
|
230 |
|
|
$bigendian = 0x55 == unpack "S", "\0U"; |
231 |
|
|
|
232 |
|
|
sub FORMAT_ULAW_BYTE (){ FORMAT_8_MULAW } |
233 |
|
|
sub FORMAT_ALAW_BYTE (){ FORMAT_8_ALAW } |
234 |
|
|
sub FORMAT_LINEAR_BYTE (){ FORMAT_8_UNSIGNED } |
235 |
|
|
sub FORMAT_LINEAR_SHORT (){ $bigendian ? FORMAT_16_LINEAR : FORMAT_16_LINEAR_LITTLE_ENDIAN } |
236 |
|
|
sub FORMAT_LINEAR_USHORT(){ $bigendian ? FORMAT_16_UNSIGNED : FORMAT_16_UNSIGNED_LITTLE_ENDIAN } |
237 |
|
|
sub FORMAT_LINEAR_LONG (){ $bigendian ? FORMAT_32_LINEAR : FORMAT_32_LINEAR_LITTLE_ENDIAN } |
238 |
|
|
sub FORMAT_LINEAR_FLOAT (){ $bigendian ? FORMAT_32_FLOAT : FORMAT_32_FLOAT_LITTLE_ENDIAN } |
239 |
|
|
sub FORMAT_LINEAR_DOUBLE(){ $bigendian ? FORMAT_32_DOUBLE : FORMAT_32_DOUBLE_LITTLE_ENDIAN } |
240 |
|
|
|
241 |
|
|
sub FILE_AU (){ FILE_NEXT } |
242 |
|
|
|
243 |
|
|
# provide some accessor methods |
244 |
|
|
for my $field (qw(rate comment filetype path format)) { |
245 |
|
|
*{"PDL::$field"} = sub { |
246 |
|
|
my $pdl = shift; |
247 |
|
|
my $hdr = $pdl->gethdr; |
248 |
|
|
if (@_) { |
249 |
|
|
$hdr->{$field} = $_[0]; |
250 |
|
|
$pdl->sethdr($hdr); |
251 |
|
|
}; |
252 |
|
|
$hdr->{$field}; |
253 |
|
|
}; |
254 |
|
|
} |
255 |
|
|
|
256 |
|
|
=head2 sound_format_name format_code |
257 |
|
|
|
258 |
|
|
Return the human-readable name of the file format with code C<format_code>. |
259 |
|
|
|
260 |
|
|
=head2 sound_type_name type_code |
261 |
|
|
|
262 |
|
|
Return the human-readable name of the sample type with code C<type_code>. |
263 |
|
|
|
264 |
|
|
=head2 describe_audio piddle |
265 |
|
|
|
266 |
|
|
Describe the audio stream contained in piddle and return it as a string. A |
267 |
|
|
fresh piddle might return: |
268 |
|
|
|
269 |
|
|
mono sound with 27411 samples |
270 |
|
|
|
271 |
|
|
Whereas a freshly loaded soundfile might yield: |
272 |
|
|
|
273 |
|
|
stereo sound with 27411 samples, original name "kongas.wav", type 2 (RIFF), |
274 |
|
|
rate 11025/s (duration 2.49s), format 7 (8-bit unsigned) |
275 |
|
|
|
276 |
|
|
=cut |
277 |
|
|
|
278 |
|
|
sub describe_audio($) { |
279 |
|
|
my $pdl = shift; |
280 |
root |
1.8 |
my ($channels, $samples) = $pdl->dims; |
281 |
root |
1.1 |
my $chan = $channels < 2 ? "mono" : |
282 |
|
|
$channels == 2 ? "stereo" : |
283 |
|
|
$channels == 4 ? "quad channel" : |
284 |
root |
1.8 |
"$channels channel"; |
285 |
root |
1.1 |
my $desc = "$chan sound with $samples samples"; |
286 |
|
|
$desc .= sprintf ", original name \"%s\"", $pdl->path if $pdl->path; |
287 |
|
|
$desc .= sprintf ", type %d (%s)", $pdl->filetype, sound_type_name($pdl->filetype) if $pdl->filetype; |
288 |
|
|
$desc .= sprintf ", rate %d/s (duration %.2fs)", $pdl->rate, $samples/$pdl->rate if $pdl->rate; |
289 |
|
|
$desc .= sprintf ", format %d (%s)", $pdl->format, sound_format_name($pdl->format) if $pdl->format; |
290 |
root |
1.8 |
$desc |
291 |
root |
1.1 |
} |
292 |
|
|
|
293 |
|
|
=head2 raudio path, [option-hash], option => value, ... |
294 |
|
|
|
295 |
root |
1.8 |
Reads audio data into the piddle. Options can be anything, most useful |
296 |
|
|
values are C<filetype>, C<rate>, C<channels> and C<format>. The returned |
297 |
|
|
piddle is represents "time" in the outer dimension, and samples in the |
298 |
root |
1.9 |
inner (i.e. scalars for mono files, 2-vectors for stereo files): |
299 |
|
|
|
300 |
|
|
[ [left0, right0], [left1, right1], [left2, right2], ...] |
301 |
root |
1.1 |
|
302 |
|
|
# read any file |
303 |
|
|
$pdl = raudio "file.wav"; |
304 |
|
|
# read a file. if it is a raw file preset values |
305 |
|
|
$pdl = raudio "file.raw", filetype => FILE_RAW, rate => 44100, channels => 2; |
306 |
|
|
|
307 |
|
|
=head2 waudio pdl, [option-hash], option => value, ... |
308 |
|
|
|
309 |
root |
1.9 |
Writes a pdl as a file. See L<raudio> for options and format. The path and |
310 |
|
|
other metadata is taken from the header, whcih cna be overwritten using |
311 |
|
|
options, e.g.: |
312 |
root |
1.1 |
|
313 |
|
|
# write a file, using the header of another piddle |
314 |
root |
1.8 |
$pdl->waudio ($orig_file->gethdr); |
315 |
root |
1.9 |
# write pdl as .au file, take rate from the header |
316 |
root |
1.8 |
$pdl->waudio (path => "piddle.au", filetype => FILE_AU, format => FORMAT_16_LINEAR; |
317 |
root |
1.1 |
|
318 |
|
|
=cut |
319 |
|
|
|
320 |
|
|
# read a sound file |
321 |
|
|
sub raudio { |
322 |
|
|
my $path = shift; |
323 |
|
|
my %hdr = (); %hdr = (%hdr, %{+shift}) if ref $_[0]; %hdr = (%hdr, @_) if @_; |
324 |
|
|
|
325 |
|
|
# for raw files this is necessary |
326 |
|
|
mus_set_raw_header_defaults $hdr{rate} || 44100, |
327 |
|
|
$hdr{channels}|| 1, |
328 |
|
|
$hdr{format} || FORMAT_LINEAR_SHORT; |
329 |
|
|
|
330 |
|
|
$hdr{path} = $path; |
331 |
|
|
$hdr{filetype} = sound_header_type $path; |
332 |
|
|
$hdr{filetype} >= 0 or barf "$path: ".audio_error_name audio_error; |
333 |
|
|
$hdr{rate} = sound_srate $path; |
334 |
|
|
$hdr{comment} = sound_comment $path if defined sound_comment $path; |
335 |
|
|
$hdr{format} = sound_data_format $path; |
336 |
|
|
my $channels = sound_chans $path; |
337 |
|
|
my $frames = sound_frames $path; |
338 |
|
|
|
339 |
|
|
my $fd = open_sound_input $path; $fd >= 0 or barf "$path: ".audio_error_name audio_error; |
340 |
|
|
|
341 |
|
|
my $pdl = zeroes long, $frames, $channels; |
342 |
|
|
$pdl->clump(-1)->read_sound ($fd, $channels, $frames) |
343 |
|
|
>= 0 or barf "$path: ".audio_error_name audio_error; |
344 |
|
|
(close_sound_input $fd) >= 0 or barf "$path: ".audio_error_name audio_error; |
345 |
|
|
$pdl = $pdl->short->xchg(0,1); |
346 |
|
|
$pdl = $pdl->clump(2) if $channels == 1; |
347 |
|
|
$pdl->sever; |
348 |
|
|
$pdl->sethdr(\%hdr); |
349 |
root |
1.8 |
$pdl |
350 |
root |
1.1 |
} |
351 |
|
|
|
352 |
|
|
sub _audio_make_plain { |
353 |
|
|
my $pdl = shift; |
354 |
|
|
if ($pdl->getndims == 1) { |
355 |
root |
1.8 |
($pdl, 1, $pdl->getdim(0)) |
356 |
root |
1.1 |
} else { |
357 |
root |
1.8 |
($pdl->xchg(0,1)->clump(-1), $pdl->dims) |
358 |
root |
1.1 |
} |
359 |
|
|
} |
360 |
|
|
|
361 |
|
|
sub waudio { |
362 |
|
|
my $pdl = shift; |
363 |
|
|
my %hdr = %{$pdl->gethdr || {}}; %hdr = (%hdr, %{+shift}) if ref $_[0]; %hdr = (%hdr, @_) if @_; |
364 |
|
|
my ($channels, $frames); |
365 |
|
|
|
366 |
|
|
$hdr{filetype} = FILE_NEXT unless exists $hdr{filetype}; |
367 |
|
|
$hdr{format} ||= FORMAT_16_LINEAR; |
368 |
|
|
$hdr{rate} ||= 44100; |
369 |
|
|
|
370 |
|
|
($pdl, $channels, $frames) = _audio_make_plain $pdl->convert(long); |
371 |
|
|
|
372 |
root |
1.8 |
1 <= $channels && $channels <= 2 |
373 |
root |
1.9 |
or croak "can only write mono or stereo (one or two channel) files, not $channels channel files"; |
374 |
root |
1.8 |
|
375 |
root |
1.1 |
my $fd = open_sound_output $hdr{path}, $hdr{rate}, $channels, $hdr{format}, $hdr{filetype}, $hdr{comment}; |
376 |
|
|
$fd >= 0 or barf "$hdr{$path}: ".audio_error_name audio_error; |
377 |
|
|
$pdl->clump(-1)->write_sound($fd, $channels, $frames) |
378 |
|
|
>= 0 or barf "$path: ".audio_error_name audio_error; |
379 |
root |
1.8 |
(close_sound_output $fd, mus_samples2bytes $hdr{format}, $frames * $channels) |
380 |
root |
1.1 |
>= 0 or barf "$hdr{$path}: ".audio_error_name audio_error; |
381 |
|
|
} |
382 |
|
|
|
383 |
|
|
=head2 cut_leading_silence pdl, level |
384 |
|
|
|
385 |
|
|
Cuts the leading silence (i.e. all samples with absolute value < level) |
386 |
|
|
and returns the resulting part. |
387 |
|
|
|
388 |
|
|
=head2 cut_trailing_silence pdl, level |
389 |
|
|
|
390 |
|
|
Cuts the trailing silence. |
391 |
|
|
|
392 |
|
|
=head2 cut_silence pdl, level |
393 |
|
|
|
394 |
|
|
Calls C<cut_leading_silence> and C<cut_trailing_silence> and returns the |
395 |
|
|
result. |
396 |
|
|
|
397 |
|
|
=cut |
398 |
|
|
|
399 |
|
|
sub cut_leading_silence { |
400 |
|
|
my $pdl = shift; |
401 |
|
|
my $level = 1*shift; |
402 |
|
|
my $skip = which (abs($pdl) > $level); |
403 |
|
|
$skip = $skip->nelem ? $skip->at(0) : 0; |
404 |
root |
1.8 |
$pdl->slice("$skip:-1") |
405 |
root |
1.1 |
} |
406 |
|
|
|
407 |
|
|
sub cut_trailing_silence { |
408 |
|
|
my $pdl = shift; |
409 |
|
|
my $level = 1*shift; |
410 |
|
|
$level = 400000; |
411 |
|
|
my $skip = which (abs($pdl) > $level); |
412 |
|
|
$skip = $skip->nelem ? $skip->at(-1) : -1; |
413 |
|
|
$skip-- if $skip > 0; |
414 |
root |
1.8 |
$pdl->slice("0:$skip") |
415 |
root |
1.1 |
} |
416 |
|
|
|
417 |
|
|
sub cut_silence { |
418 |
|
|
$_[0]->cut_leading_silence($_[1]) |
419 |
root |
1.8 |
->cut_trailing_silence($_[1]) |
420 |
root |
1.1 |
} |
421 |
|
|
|
422 |
|
|
# have we been a bad boy? |
423 |
|
|
|
424 |
|
|
for (@METHODS) { |
425 |
|
|
*{"PDL::$_"} = \&$_; |
426 |
|
|
} |
427 |
|
|
EOD |
428 |
|
|
|
429 |
|
|
for my $d (qw(read write)) { |
430 |
|
|
pp_def($d.'_sound', |
431 |
|
|
Pars => ($d eq "read" ? '[o]' : '') . 'sample(n)', |
432 |
|
|
OtherPars => "int fd; int chans; int frames", |
433 |
|
|
GenericTypes => [L], |
434 |
|
|
PMCode => '', |
435 |
|
|
Doc => undef, |
436 |
|
|
Code => q@ |
437 |
|
|
int chans = $COMP(chans); |
438 |
|
|
int frames = $COMP(frames); |
439 |
|
|
int **bufs = malloc (chans * sizeof (int *)); |
440 |
|
|
int i; |
441 |
|
|
dSP; |
442 |
|
|
|
443 |
|
|
/* hmm... */ |
444 |
|
|
assert (sizeof (PDL_long) == sizeof (int)); |
445 |
|
|
|
446 |
|
|
for (i = 0; i < chans; i++) |
447 |
|
|
bufs[i] = (int *)$P(sample) + i * frames; |
448 |
|
|
|
449 |
|
|
XPUSHs (sv_2mortal (newSViv (@.$d.q@_sound ($COMP(fd), 0, frames - 1, chans, bufs)))); |
450 |
|
|
|
451 |
|
|
free (bufs); |
452 |
|
|
@, |
453 |
|
|
); |
454 |
|
|
} |
455 |
|
|
|
456 |
|
|
pp_addhdr <<'EOH'; |
457 |
|
|
#include "sndlib/sndlib.h" |
458 |
|
|
#include "xlib.h" |
459 |
|
|
|
460 |
|
|
static void |
461 |
|
|
mus_barfer (int err_type, char *err_msg) |
462 |
|
|
{ |
463 |
|
|
barf ("%s [MUSERR]", err_msg); |
464 |
|
|
} |
465 |
|
|
|
466 |
|
|
EOH |
467 |
|
|
|
468 |
|
|
pp_addpm <<'EOP'; |
469 |
|
|
|
470 |
|
|
=head2 playaudio pdl, [option-hash], option => value ... |
471 |
|
|
|
472 |
|
|
Play the piddle as an audio file. Options can be supplied either through |
473 |
|
|
the option hash (a hash-reference), through the pdl header or the options: |
474 |
|
|
|
475 |
|
|
# play a piddle that has a valid header (e.g. from raudio) |
476 |
|
|
$pdl->playaudio; |
477 |
|
|
# play it with a different samplerate |
478 |
|
|
$pdl->playaudio(rate => 22050); |
479 |
|
|
|
480 |
|
|
=cut |
481 |
|
|
|
482 |
|
|
EOP |
483 |
|
|
|
484 |
|
|
pp_def('playaudio', |
485 |
|
|
Pars => 'sample(n)', |
486 |
|
|
OtherPars => "int srate; int chans; int format; int dev; int bufsize", |
487 |
|
|
GenericTypes => [B,U,S,L], |
488 |
|
|
Doc => undef, |
489 |
|
|
PMCode => q! |
490 |
|
|
sub PDL::playaudio { |
491 |
|
|
my $pdl = shift; |
492 |
|
|
my %hdr = %{$pdl->gethdr || {}}; %hdr = (%hdr, %{+shift}) if ref $_[0]; %hdr = (%hdr, @_) if @_; |
493 |
|
|
my $chans; |
494 |
|
|
$hdr{format} ||= FORMAT_LINEAR_SHORT; |
495 |
|
|
($pdl, $chans) = _audio_make_plain $pdl; |
496 |
|
|
&PDL::_playaudio_int($pdl, |
497 |
|
|
$hdr{rate} || 44100, |
498 |
|
|
$chans, |
499 |
|
|
$hdr{format}, |
500 |
|
|
$hdr{device} || DEV_DEFAULT, |
501 |
|
|
$hdr{bufsize} || 32768); |
502 |
|
|
} |
503 |
|
|
!, |
504 |
|
|
Code => q@ |
505 |
|
|
int ns; |
506 |
|
|
int line; |
507 |
|
|
|
508 |
|
|
if (sizeof ($GENERIC ()) != mus_format2bytes ($COMP(format))) |
509 |
|
|
barf ("pdl datatype and selected audio format do not match"); |
510 |
|
|
|
511 |
|
|
if (initialize_audio () < 0) |
512 |
|
|
barf ("playaudio: %s", audio_error_name (audio_error ())); |
513 |
|
|
|
514 |
|
|
if ((line = open_audio_output ($COMP(dev), $COMP(srate), $COMP(chans), $COMP(format), $COMP(bufsize))) < 0) |
515 |
|
|
barf ("playaudio: %s", audio_error_name (audio_error ())); |
516 |
|
|
|
517 |
|
|
ns = $SIZE(n); |
518 |
|
|
threadloop %{ |
519 |
|
|
if (write_audio (line, (void *)$P(sample), ns * sizeof ($GENERIC()))) |
520 |
|
|
barf ("playaudio: %s", audio_error_name (audio_error ())); |
521 |
|
|
%} |
522 |
|
|
if (close_audio (line) < 0) |
523 |
|
|
barf ("playaudio: %s", audio_error_name (audio_error ())); |
524 |
|
|
@, |
525 |
|
|
); |
526 |
|
|
|
527 |
|
|
pp_def 'ulaw2linear', |
528 |
|
|
Pars => 'byte u(n); short [o] s(n)', |
529 |
|
|
GenericTypes => [B], |
530 |
|
|
Doc => 'conversion from (m)u-law into signed, linear, 16 bit samples (rather slow)', |
531 |
|
|
Code => q! loop(n) %{ $s() = st_ulaw_to_linear ($u()); %} !, |
532 |
|
|
; |
533 |
|
|
|
534 |
|
|
pp_def 'linear2ulaw', |
535 |
|
|
Pars => 'short s(n); byte [o] u(n)', |
536 |
|
|
GenericTypes => [S], |
537 |
|
|
Doc => 'conversion from signed, linear, 16 bit samples into (m)u-law (rather slow)', |
538 |
|
|
Code => q! loop(n) %{ $u() = st_linear_to_ulaw ($s()); %} !, |
539 |
|
|
; |
540 |
|
|
|
541 |
|
|
pp_def 'alaw2linear', |
542 |
|
|
Pars => 'byte u(n); short [o] s(n)', |
543 |
|
|
GenericTypes => [B], |
544 |
|
|
Doc => 'conversion from A-law into signed, linear, 16 bit samples (rather slow)', |
545 |
|
|
Code => q! loop(n) %{ $s() = st_Alaw_to_linear ($u()); %} !, |
546 |
|
|
; |
547 |
|
|
|
548 |
|
|
pp_def 'linear2alaw', |
549 |
|
|
Pars => 'short s(n); byte [o] u(n)', |
550 |
|
|
GenericTypes => [S], |
551 |
|
|
Doc => 'conversion from signed, linear, 16 bit samples into A-law (rather slow)', |
552 |
|
|
Code => q! loop(n) %{ $u() = st_linear_to_Alaw ($s()); %} !, |
553 |
|
|
; |
554 |
|
|
|
555 |
|
|
pp_addpm <<'EOP'; |
556 |
|
|
|
557 |
|
|
=head2 gen_oscil duration*, freq/, phase-mod, [fm-mod/] |
558 |
|
|
|
559 |
|
|
=head2 gen_sawtooth duration*, freq/, phase-mod, [fm-mod/] |
560 |
|
|
|
561 |
|
|
=head2 gen_square duration*, freq/, phase-mod, duty, [fm-mod/] |
562 |
|
|
|
563 |
|
|
=head2 gen_triangle duration*, freq/, phase-mod, [fm-mod/] |
564 |
|
|
|
565 |
|
|
=head2 gen_pulse_train duration*, freq/, phase-mod, [fm-mod/] |
566 |
|
|
|
567 |
|
|
=head2 gen_rand duration*, freq/ |
568 |
|
|
|
569 |
|
|
=head2 gen_rand_1f duration* |
570 |
|
|
|
571 |
|
|
All of these functions generate appropriate waveforms with frequency |
572 |
|
|
C<freq> (cycles/sample) and phase C<phase> (0..1). |
573 |
|
|
|
574 |
|
|
The C<duration> might be either a piddle (which gives the form of the |
575 |
|
|
output) or the number of samples to generate. |
576 |
|
|
|
577 |
|
|
The output samples are between -1 and +1 (i.e. "-1 <= s <= +1"). |
578 |
|
|
|
579 |
|
|
The C<duty> parameter of the square generator influences the duty cycle of |
580 |
|
|
the signal. Zero means 50%-50%, 0.5 means 75% on, 25% off, -0.8 means 10% |
581 |
|
|
on, 90% off etc... Of course, the C<duty> parameter might also be a vector |
582 |
|
|
of size C<duration>. |
583 |
|
|
|
584 |
|
|
=cut |
585 |
|
|
|
586 |
|
|
# return the time offset for duration, freq, phase, fm_mod |
587 |
|
|
sub _dur2time($$;$$) { |
588 |
|
|
my ($dur, $freq, $phase, $fm_mod) = @_; |
589 |
|
|
zeroes($dur)->xvals*$freq+cumusumover($fm_mod)+$phase |
590 |
|
|
} |
591 |
|
|
|
592 |
|
|
sub gen_oscil($$;$$) { |
593 |
|
|
my ($dur, $freq, $phase, $fm_mod) = @_; |
594 |
|
|
sin M_2PI*_dur2time($dur,$freq,$phase,$fm_mod); |
595 |
|
|
} |
596 |
|
|
|
597 |
|
|
sub gen_sawtooth($$;$$) { |
598 |
|
|
my ($dur, $freq, $phase, $fm_mod) = @_; |
599 |
|
|
# maybe better use abs? |
600 |
|
|
$freq = $fm_mod * $freq if defined $fm_mod; |
601 |
|
|
(_dur2time($dur,$freq,$phase,$fm_mod) % 1 - 0.5) * 2; |
602 |
|
|
} |
603 |
|
|
|
604 |
|
|
sub gen_square($$;$$$) { |
605 |
|
|
my ($dur, $freq, $phase, $duty, $fm_mod) = @_; |
606 |
|
|
(gen_sawtooth($dur,$freq,$phase,$fm_mod) < $duty) * 2 - 1; |
607 |
|
|
} |
608 |
|
|
|
609 |
|
|
sub gen_triangle($$;$$) { |
610 |
|
|
my ($dur, $freq, $phase, $fm_mod) = @_; |
611 |
|
|
my $st = gen_sawtooth($dur,$freq,$phase,$fm_mod); |
612 |
|
|
$st * (($st < 0) * 2 - 1); |
613 |
|
|
} |
614 |
|
|
|
615 |
|
|
sub gen_pulse_train($$;$$) { |
616 |
|
|
my ($dur, $freq, $phase, $fm_mod) = @_; |
617 |
|
|
my $st = gen_sawtooth($dur,$freq,$phase,$fm_mod); |
618 |
|
|
$st < $st->rotate(1); |
619 |
|
|
} |
620 |
|
|
|
621 |
|
|
sub gen_rand($$) { |
622 |
|
|
my ($dur, $freq) = @_; |
623 |
|
|
my $z = zeroes $dur; |
624 |
|
|
my $r = (random zeroes $freq * $z->getdim(0) + 1) * 2 - 1; |
625 |
|
|
$r->index($freq * $z->xvals)->sever; |
626 |
|
|
} |
627 |
|
|
|
628 |
|
|
sub gen_rand_1f($) { |
629 |
|
|
my ($dur) = @_; |
630 |
|
|
$dur = zeroes $dur; |
631 |
|
|
_gen_noise_1f($dur); |
632 |
|
|
$dur; |
633 |
|
|
} |
634 |
|
|
|
635 |
|
|
=head2 gen_env duration*, xvals, yvals, [base] |
636 |
|
|
|
637 |
|
|
Generates an interpolated envelope between the points given by xvals and |
638 |
|
|
yvals. When base == 1 (the default) then the values will be linearly |
639 |
|
|
interpolated, otherwise they follow an exponential curve that is bend |
640 |
|
|
inwards (base < 1) or outwards (base > 1). |
641 |
|
|
|
642 |
|
|
# generate a linear envelope with attack in the first 10% |
643 |
|
|
gen_env 5000, [0 1 2 9 10], [0 1 0.6 0.6 0]; |
644 |
|
|
|
645 |
|
|
=cut |
646 |
|
|
|
647 |
|
|
sub gen_env($;@) { |
648 |
|
|
my ($dur, $x, $y, $base) = @_; |
649 |
|
|
my $pdl = zeroes $dur; |
650 |
|
|
$base ||= 1; |
651 |
|
|
$base == 1 or barf "base values != 1 are not yet implemented\n"; |
652 |
|
|
my @x = ref $x eq "ARRAY" ? @$x : $x->list; |
653 |
|
|
my @y = ref $y eq "ARRAY" ? @$y : $y->list; |
654 |
|
|
my $f = $pdl->getdim(0) / $x[-1]; |
655 |
|
|
@x == @y or barf "gen_env: x and y lists have different sizes"; |
656 |
|
|
my $x0 = 0; |
657 |
|
|
my $y0 = 0; |
658 |
|
|
for (0..$#x) { |
659 |
|
|
my ($x,$y) = (int($x[$_]*$f),$y[$_]); |
660 |
|
|
if ($x > $x0) { |
661 |
|
|
my $a = $pdl->slice("$x0:".($x-1)); |
662 |
|
|
if ($y0 != $y) { |
663 |
|
|
$a .= $a->xlinvals($y0,$y); |
664 |
|
|
} else { |
665 |
|
|
$a .= $y; |
666 |
|
|
} |
667 |
|
|
} |
668 |
|
|
($x0,$y0)=($x,$y); |
669 |
|
|
} |
670 |
|
|
$pdl; |
671 |
|
|
} |
672 |
|
|
|
673 |
|
|
=head2 gen_adsr duration*, sustain-level, attack-time, decay-time, sustain-time, release-time |
674 |
|
|
|
675 |
|
|
Simple ADSR envelope generator. The C<sustain-level> is the amplitude (0 |
676 |
|
|
to 1) of the sustain level. The other for parameters give the relative |
677 |
|
|
interval times, in any unit you like, only their relative ratios |
678 |
|
|
are important. Any of these times might be zero, in which case the |
679 |
|
|
corresponding part is omitted from the envelope. |
680 |
|
|
|
681 |
|
|
=cut |
682 |
|
|
|
683 |
|
|
sub gen_adsr($$$$$$) { |
684 |
|
|
my ($dur, $amp, $a, $d, $s, $r) = @_; |
685 |
|
|
my (@x, @y); |
686 |
|
|
my $t = 0; |
687 |
|
|
push @x, 0; push @y, 0; |
688 |
|
|
if ($a > 0) { push @x, $t += $a; push @y, 1 } else { $y[-1] = 1 } |
689 |
|
|
if ($d > 0) { push @x, $t += $d; push @y, $amp } else { $y[-1] = $amp } |
690 |
|
|
if ($s > 0) { push @x, $t += $s; push @y, $amp } else { $y[-1] = $amp } |
691 |
|
|
if ($r > 0) { push @x, $t += $r; push @y, 0 } else { $y[-1] = 0 } |
692 |
|
|
zeroes($dur)->xlinvals(0,$t)->linear_interpolate(pdl(@x), pdl(@y)); |
693 |
|
|
} |
694 |
|
|
|
695 |
|
|
=head2 gen_asymmetric_fm duration*, freq/, phase, [r , [ratio]] |
696 |
|
|
|
697 |
|
|
C<gen_asymmetric_fm> provides a way around the symmetric spectra normally |
698 |
|
|
produced by FM. See Palamin and Palamin, "A Method of Generating and |
699 |
|
|
Controlling Asymmetrical Spectra" JAES vol 36, no 9, Sept 88, p671-685. |
700 |
|
|
|
701 |
|
|
=cut |
702 |
|
|
|
703 |
|
|
sub gen_asymmetric_fm($$;$$$) { |
704 |
|
|
my ($dur, $freq, $phase, $r, $ratio) = @_; |
705 |
|
|
$r ||= 1; |
706 |
|
|
$ratio ||= 1; |
707 |
|
|
my $cosr = 0.5 * ($r - 1/$r); |
708 |
|
|
my $sinr = 0.5 * ($r + 1/$r); |
709 |
|
|
$phase = zeroes($dur)->xvals*$freq + $phase; |
710 |
|
|
my $mth = $phase * $ratio; |
711 |
|
|
exp ($cosr * cos $mth) * sin ($phase + $sinr * sin $mth); |
712 |
|
|
} |
713 |
|
|
|
714 |
|
|
=head2 gen_sum_of_cosines duration*, freq/, phase, ncosines, [fm_mod/] |
715 |
|
|
|
716 |
|
|
Generates a sum of C<n> cosines C<(1 + 2(cos(x) + cos(2x) + ... cos(nx)) = |
717 |
|
|
sin((n+.5)x) / sin(x/2))>. Other arguments are similar to to C<gen_oscil>. |
718 |
|
|
|
719 |
|
|
=head2 gen_sine_summation duration*, freq/, phase, [nsines, [a, [b_ratio, [fm_mod/]]]] |
720 |
|
|
|
721 |
|
|
C<gen_sine_summation> provides a kind of additive synthesis. See |
722 |
|
|
J.A.Moorer, "Signal Processing Aspects of Computer Music" and "The |
723 |
|
|
Synthesis of Complex Audio Spectra by means of Discrete Summation |
724 |
|
|
Formulae" (Stan-M-5). The basic idea is very similar to that used in |
725 |
|
|
gen_sum_of_cosines generator. |
726 |
|
|
|
727 |
|
|
The default value for C<nsines> is 1 (but zero is a valid value), for C<a> |
728 |
|
|
is 0.5 and for C<b_ratio> is 1. |
729 |
|
|
|
730 |
|
|
(btw, either my formula is broken or the output indeed does not lie |
731 |
|
|
between -1 and +1, but rather -5 .. +5). |
732 |
|
|
|
733 |
|
|
=cut |
734 |
|
|
|
735 |
|
|
sub gen_sum_of_cosines { |
736 |
|
|
my ($dur, $freq, $phase, $cosines, $fm_mod) = @_; |
737 |
|
|
$cosines ||= 1; |
738 |
|
|
$dur = M_2PI*_dur2time($dur,$freq,$phase+1e-6,$fm_mod); |
739 |
|
|
sin ($dur * ($cosines+0.5)) / (sin ($dur * 0.5) * (1+2*$cosines)); |
740 |
|
|
} |
741 |
|
|
|
742 |
|
|
sub gen_sine_summation { |
743 |
|
|
my ($dur, $freq, $phase, $n, $a, $ratio, $fm_mod) = @_; |
744 |
|
|
$n = 1 unless defined $n; |
745 |
|
|
$a ||= 0.5; |
746 |
|
|
$ratio ||= 1; |
747 |
|
|
|
748 |
|
|
$dur = M_2PI*_dur2time($dur,$freq,$phase+1e-6,$fm_mod); |
749 |
|
|
my $an = $n ? $a**($n+1) : 0; |
750 |
|
|
my $B = $dur * $ratio; |
751 |
|
|
(sin ($dur) |
752 |
|
|
- $a * sin ($dur - $B) |
753 |
|
|
- $an * (sin ($dur + $B*($n+1)) |
754 |
|
|
- $a * sin ($dur + $B*$n))) |
755 |
|
|
/ (1 + $a*$a - (2*$a*cos($B))); |
756 |
|
|
} |
757 |
|
|
|
758 |
|
|
=head2 gen_from_table duration*, frequency/, table, [phase], [fm_mod/] |
759 |
|
|
|
760 |
|
|
C<gen_from_table> generates a waveform by repeating a waveform given |
761 |
|
|
in C<table>, linearly interpolating between successive points of the |
762 |
|
|
C<waveform>. |
763 |
|
|
|
764 |
|
|
=head2 partials2waveshape size*, partials, amplitudes, [phase], [fm_mod/] |
765 |
|
|
|
766 |
root |
1.6 |
Take a (perl or pdl) list of (integer) C<partials> and a list of |
767 |
root |
1.1 |
C<amplitudes> and generate a single wave shape that results by adding |
768 |
|
|
these partial sines. |
769 |
|
|
|
770 |
|
|
This could (and should) be used by the C<gen_from_table> generator. |
771 |
|
|
|
772 |
|
|
=head2 gen_from_partials duration*, frequency/, partials, amplitudes, [phase], [fm_mod/] |
773 |
|
|
|
774 |
root |
1.6 |
Take a (perl list or pdl) list of (possibly noninteger) C<partials> and a |
775 |
root |
1.1 |
list of C<amplitudes> and generate the waveform resulting by summing up |
776 |
|
|
all these partial sines. |
777 |
|
|
|
778 |
|
|
=cut |
779 |
|
|
|
780 |
|
|
sub gen_from_table { |
781 |
|
|
my ($dur, $freq, $table, $phase, $fm_mod) = @_; |
782 |
|
|
my $r = _dur2time($dur,$freq,$phase,$fm_mod); |
783 |
|
|
my $tx = ($table->dims)[0]; |
784 |
|
|
$table = cat($table, $table->at(0))->clump(-1); |
785 |
|
|
$r %= 1; $r += 1; $r %= 1; $r *= $tx; |
786 |
|
|
$r->linear_interpolate($table->xvals, $table); |
787 |
|
|
} |
788 |
|
|
|
789 |
|
|
sub partials2waveshape { |
790 |
|
|
my $dur = zeroes $_[0]; |
791 |
|
|
my @idx = ref $_[1] eq "ARRAY" ? @{$_[1]} : $_[1]->list; |
792 |
|
|
my @amp = ref $_[2] eq "ARRAY" ? @{$_[2]} : $_[2]->list; |
793 |
|
|
|
794 |
|
|
my $angle = $dur->xvals * (M_2PI/$dur->getdim(0)); |
795 |
|
|
|
796 |
|
|
for (0..$#idx) { |
797 |
|
|
$dur += $amp[$_] * sin $idx[$_]*$angle; |
798 |
|
|
} |
799 |
|
|
$dur / max abs $dur; |
800 |
|
|
} |
801 |
|
|
|
802 |
|
|
sub gen_from_partials { |
803 |
|
|
my ($dur, $freq, $ti, $ta, $phase, $fm_mod) = @_; |
804 |
|
|
my @idx = ref $ti eq "ARRAY" ? @{$ti} : $ti->list; |
805 |
|
|
my @amp = ref $ta eq "ARRAY" ? @{$ta} : $ta->list; |
806 |
|
|
|
807 |
|
|
$dur = zeroes $dur; |
808 |
|
|
|
809 |
|
|
my $angle = _dur2time($dur,$freq,$phase,$fm_mod); |
810 |
|
|
|
811 |
|
|
for (0..$#idx) { |
812 |
|
|
$dur += $amp[$_] * fast_sin($idx[$_]*$angle); |
813 |
|
|
} |
814 |
|
|
$dur / max abs $dur; |
815 |
|
|
} |
816 |
|
|
|
817 |
|
|
EOP |
818 |
|
|
|
819 |
|
|
pp_def '_gen_noise_1f', |
820 |
|
|
Pars => '[o]f(n)', |
821 |
|
|
GenericTypes => [F,D], |
822 |
|
|
Doc => undef, |
823 |
|
|
PMCode => undef, |
824 |
|
|
Code => q^ |
825 |
|
|
static unsigned long ctz[64] = { |
826 |
|
|
6, 0, 1, 0, 2, 0, 1, 0, |
827 |
|
|
3, 0, 1, 0, 2, 0, 1, 0, |
828 |
|
|
4, 0, 1, 0, 2, 0, 1, 0, |
829 |
|
|
3, 0, 1, 0, 2, 0, 1, 0, |
830 |
|
|
5, 0, 1, 0, 2, 0, 1, 0, |
831 |
|
|
3, 0, 1, 0, 2, 0, 1, 0, |
832 |
|
|
4, 0, 1, 0, 2, 0, 1, 0, |
833 |
|
|
3, 0, 1, 0, 2, 0, 1, 0, |
834 |
|
|
}; |
835 |
|
|
static unsigned int dice[7]; |
836 |
|
|
static unsigned int total = 0; |
837 |
|
|
unsigned int n, prevrand, newrand, seed, k; |
838 |
|
|
|
839 |
|
|
threadloop %{ |
840 |
|
|
for (n = 0; n < $SIZE(n); n++) |
841 |
|
|
{ |
842 |
|
|
k = ctz [n & 63]; |
843 |
|
|
|
844 |
|
|
prevrand = dice[k]; |
845 |
|
|
seed = 1664525 * seed + 1013904223; |
846 |
|
|
newrand = (seed >> 3); |
847 |
|
|
dice[k] = newrand; |
848 |
|
|
|
849 |
|
|
total += newrand - prevrand; |
850 |
|
|
|
851 |
|
|
seed = 1103515245 * seed + 12345; |
852 |
|
|
newrand = (seed >> 3); |
853 |
|
|
|
854 |
|
|
$f(n=>n) = ($GENERIC())(total + newrand) * (1./(3<<29)) - 1; |
855 |
|
|
} |
856 |
|
|
%} |
857 |
|
|
^ |
858 |
|
|
; |
859 |
|
|
|
860 |
|
|
pp_def 'filter_one_zero', |
861 |
|
|
Pars => 'in(n); [o] out(n)', |
862 |
|
|
OtherPars => 'double a0; double a1', |
863 |
|
|
GenericTypes => [F,D], |
864 |
|
|
Doc => 'apply a one zero filter, y(n) = a0 x(n) + a1 x(n-1)', |
865 |
|
|
Code => q! |
866 |
|
|
double a0 = $COMP(a0); |
867 |
|
|
double a1 = $COMP(a1); |
868 |
|
|
double x1 = 0.; |
869 |
|
|
|
870 |
|
|
threadloop %{ loop(n) %{ |
871 |
|
|
$out() = a0 * $in() + a1 * x1; |
872 |
|
|
x1 = $in(); |
873 |
|
|
%} %} |
874 |
|
|
! |
875 |
|
|
; |
876 |
|
|
|
877 |
|
|
pp_def 'filter_one_pole', |
878 |
|
|
Pars => 'in(n); [o] out(n)', |
879 |
|
|
OtherPars => 'double a0; double b1', |
880 |
|
|
GenericTypes => [F,D], |
881 |
|
|
Doc => 'apply a one pole filter, y(n) = a0 x(n) - b1 y(n-1)', |
882 |
|
|
Code => q! |
883 |
|
|
double a0 = $COMP(a0); |
884 |
|
|
double b1 = $COMP(b1); |
885 |
|
|
double y2 = 0.; |
886 |
|
|
|
887 |
|
|
threadloop %{ loop(n) %{ |
888 |
|
|
$out() = y2 = a0 * $in() - b1 * y2; |
889 |
|
|
%} %} |
890 |
|
|
! |
891 |
|
|
; |
892 |
|
|
|
893 |
|
|
pp_def 'filter_two_zero', |
894 |
|
|
Pars => 'in(n); [o] out(n)', |
895 |
|
|
OtherPars => 'double a0; double a1; double a2', |
896 |
|
|
GenericTypes => [F,D], |
897 |
|
|
Doc => 'apply a two zero filter, y(n) = a0 x(n) + a1 x(n-1) + a2 x(n-2)', |
898 |
|
|
Code => q! |
899 |
|
|
double a0 = $COMP(a0); |
900 |
|
|
double a1 = $COMP(a1); |
901 |
|
|
double a2 = $COMP(a2); |
902 |
|
|
double x1 = 0.; |
903 |
|
|
double x2 = 0.; |
904 |
|
|
|
905 |
|
|
threadloop %{ loop(n) %{ |
906 |
|
|
$out() = a0 * $in() + a1 * x1 + a2 * x2; |
907 |
|
|
x2 = x1; x1 = $in(); |
908 |
|
|
%} %} |
909 |
|
|
! |
910 |
|
|
; |
911 |
|
|
|
912 |
|
|
pp_def 'filter_two_pole', |
913 |
|
|
Pars => 'in(n); [o] out(n)', |
914 |
|
|
OtherPars => 'double a0; double b1; double b2', |
915 |
|
|
GenericTypes => [F,D], |
916 |
|
|
Doc => 'apply a two pole filter, y(n) = a0 x(n) - b1 y(n-1) - b2 y(n-2)', |
917 |
|
|
Code => q! |
918 |
|
|
double a0 = $COMP(a0); |
919 |
|
|
double b1 = $COMP(b1); |
920 |
|
|
double b2 = $COMP(b2); |
921 |
|
|
double y1 = 0.; |
922 |
|
|
double y2 = 0.; |
923 |
|
|
|
924 |
|
|
threadloop %{ loop(n) %{ |
925 |
|
|
$out() = a0 * $in() - b1 * y1 - b2 * y2; |
926 |
|
|
y2 = y1; y1 = $out(); |
927 |
|
|
%} %} |
928 |
|
|
! |
929 |
|
|
; |
930 |
|
|
|
931 |
|
|
pp_def 'filter_formant', |
932 |
|
|
Pars => 'in(n); [o] out(n)', |
933 |
|
|
OtherPars => 'double radius; double frequency; double gain', |
934 |
|
|
GenericTypes => [F,D], |
935 |
|
|
Doc => 'apply a formant filter, y(n) = x(n) - r*x(n-2) + 2*r*cos(2*pi*frequency)*y(n-1) - r*r*y(n-2). A good value for C<gain> is 1.', |
936 |
|
|
Code => q! |
937 |
|
|
threadloop %{ |
938 |
|
|
double radius = 1. - $COMP(radius); |
939 |
|
|
double a0 = $COMP(gain) * sin ($COMP(frequency) * M_2PI) * (1. - (radius * radius)); |
940 |
|
|
double a2 = -radius; |
941 |
|
|
double b1 = -2. * radius * cos ($COMP(frequency) * M_2PI); |
942 |
|
|
double b2 = radius * radius; |
943 |
|
|
double x1 = 0.; |
944 |
|
|
double x2 = 0.; |
945 |
|
|
double y1 = 0.; |
946 |
|
|
double y2 = 0.; |
947 |
|
|
|
948 |
|
|
loop(n) %{ |
949 |
|
|
double inval = a0 * $in(); |
950 |
|
|
$out() = inval + a2 * x2 - b1 * y1 - b2 * y2; |
951 |
|
|
y2 = y1; y1 = $out(); |
952 |
|
|
x2 = x1; x1 = inval; |
953 |
|
|
%} |
954 |
|
|
%} |
955 |
|
|
! |
956 |
|
|
; |
957 |
|
|
|
958 |
|
|
pp_addpm <<'EOP'; |
959 |
|
|
|
960 |
|
|
=head2 filter_ppolar pdl, radius/, frequency/ |
961 |
|
|
|
962 |
|
|
apply a two pole filter (given in polar form). The filter has two poles, |
963 |
|
|
one at (radius,frequency), the other at (radius,-frequency). Radius |
964 |
|
|
is between 0 and 1 (but less than 1), and frequency is between 0 and |
965 |
|
|
0.5. This is the standard resonator form with poles specified by the |
966 |
|
|
polar coordinates of one pole. |
967 |
|
|
|
968 |
|
|
=head2 filter_zpolar pdl, radius/, frequency/ |
969 |
|
|
|
970 |
|
|
apply a two zero filter (given in polar form). See C<filter_ppolar>. |
971 |
|
|
|
972 |
|
|
=cut |
973 |
|
|
|
974 |
|
|
sub filter_ppolar { |
975 |
|
|
my ($pdl, $radius, $freq) = @_; |
976 |
|
|
$pdl->filter_two_pole(1, -2*$radius*cos($freq*M_2PI), $radius**2); |
977 |
|
|
} |
978 |
|
|
|
979 |
|
|
sub filter_zpolar { |
980 |
|
|
my ($pdl, $radius, $freq) = @_; |
981 |
|
|
$pdl->filter_two_zero(1, -2*$radius*cos($freq*M_2PI), $radius**2); |
982 |
|
|
} |
983 |
|
|
|
984 |
|
|
=head2 partials2polynomial partials, [kind] |
985 |
|
|
|
986 |
|
|
C<partials2polynomial> takes a list of harmonic amplitudes and returns a |
987 |
|
|
list of Chebychev polynomial coefficients. The argument C<kind> determines |
988 |
|
|
which kind of Chebychev polynomial we are interested in, 1st kind or 2nd |
989 |
|
|
kind. (default is 1). |
990 |
|
|
|
991 |
|
|
=cut |
992 |
|
|
|
993 |
|
|
sub partials2polynomial { |
994 |
|
|
my ($partials, $kind) = @_; |
995 |
|
|
$kind = 1 unless defined $kind; |
996 |
|
|
my $t0 = zeroes($partials->getdim(0)+1); |
997 |
|
|
my $t1 = $t0->copy; |
998 |
|
|
my $cc1 = $t0->copy; |
999 |
|
|
|
1000 |
|
|
$t0->set(0,1); |
1001 |
|
|
$t1->set(1,$kind); |
1002 |
|
|
|
1003 |
|
|
for my $amp ($partials->list) { |
1004 |
|
|
$cc1 += $amp * $t1; |
1005 |
|
|
$tn = 2 * $t1->rshift(1) - $t0; |
1006 |
|
|
$t0 = $t1; |
1007 |
|
|
$t1 = $tn; |
1008 |
|
|
} |
1009 |
|
|
$cc1; |
1010 |
|
|
} |
1011 |
|
|
|
1012 |
|
|
=head2 ring_modulate in1, in2 |
1013 |
|
|
|
1014 |
|
|
ring modulates in1 with in2 (this is just a multiply). |
1015 |
|
|
|
1016 |
|
|
=head2 amplitude_modulate am_carrier, in1, in2 |
1017 |
|
|
|
1018 |
|
|
amplitude modulates am_carrier and in2 with in1 (this calculates in1 * (am_carrier + in2)). |
1019 |
|
|
|
1020 |
|
|
=cut |
1021 |
|
|
|
1022 |
|
|
sub ring_modulate { |
1023 |
|
|
$_[0] * $_[1]; |
1024 |
|
|
} |
1025 |
|
|
|
1026 |
|
|
sub amplitude_modulate { |
1027 |
|
|
$_[1] * ($_[0] + $_[2]); |
1028 |
|
|
} |
1029 |
|
|
|
1030 |
|
|
EOP |
1031 |
|
|
|
1032 |
|
|
pp_def 'filter_sir', |
1033 |
|
|
Pars => 'x(n); a(an); b(bn); [o]y(n)', |
1034 |
|
|
GenericTypes => [L,F,D], |
1035 |
|
|
Doc => <<'EOD', |
1036 |
|
|
Generic (short delay) impulse response filter. C<x> is the input signal |
1037 |
|
|
(which is supposed to be zero for negative indices). C<a> contains the |
1038 |
|
|
input (x) coefficients (a0, a1, .. an), whereas C<b> contains the output |
1039 |
|
|
(y) coefficients (b0, b1, ... bn), i.e.: |
1040 |
|
|
|
1041 |
|
|
y(n) = a0 x(n) - b1 y(n-1) + a1 x(n-1) - b2 y(n-2) + a2 x(n-2) - b3 ... |
1042 |
|
|
|
1043 |
|
|
This can be used to generate fir and iir filters of any length, or even |
1044 |
|
|
more complicated constructs. |
1045 |
|
|
|
1046 |
|
|
C<b0> (then first element of C<b>) is being ignored currently AND SHOULD |
1047 |
|
|
BE SPECIFIED AS ONE FOR FUTURE COMPATIBILITY |
1048 |
|
|
EOD |
1049 |
|
|
Code => q^ |
1050 |
|
|
int an = $SIZE(an); |
1051 |
|
|
int bn = $SIZE(bn); |
1052 |
|
|
int i, n, I; |
1053 |
|
|
$GENERIC() v; |
1054 |
|
|
|
1055 |
|
|
threadloop %{ |
1056 |
|
|
for (n = 0; n < $SIZE(n); n++) |
1057 |
|
|
{ |
1058 |
|
|
v = 0; |
1059 |
|
|
|
1060 |
|
|
/* apply loop splitting manually later! */ |
1061 |
|
|
|
1062 |
|
|
for (i = 0; i < an; i++) |
1063 |
|
|
{ |
1064 |
|
|
I = n-i; |
1065 |
|
|
if (I >= 0) |
1066 |
|
|
v += $a(an=>i) * $x(n=>I); |
1067 |
|
|
} |
1068 |
|
|
|
1069 |
|
|
for (i = 0; i < bn; i++) |
1070 |
|
|
{ |
1071 |
|
|
I = n-i; |
1072 |
|
|
if (I >= 0) |
1073 |
|
|
v -= $b(bn=>i) * $y(n=>I); |
1074 |
|
|
} |
1075 |
|
|
|
1076 |
|
|
$y(n=>n) = v; |
1077 |
|
|
} |
1078 |
|
|
%} |
1079 |
|
|
^ |
1080 |
|
|
; |
1081 |
|
|
|
1082 |
|
|
pp_def 'filter_lir', |
1083 |
|
|
Pars => 'x(n); int a_x(an); a_y(an); int b_x(bn); b_y(bn); [o]y(n)', |
1084 |
|
|
GenericTypes => [L,F,D], |
1085 |
|
|
Doc => <<'EOD', |
1086 |
|
|
|
1087 |
|
|
Generic (long delay) impulse response filter. The difference to |
1088 |
|
|
C<filter_sir> is that the filter coefficients need not be consecutive, |
1089 |
|
|
but instead their indices are given by the C<a_x> and C<b_x> (integer) |
1090 |
|
|
vectors, while the corresponding coefficients are in C<a_y> and |
1091 |
|
|
C<b_y>. (All C<a_x> must be >= 0, while all the C<b_x> must be >= 1, as |
1092 |
|
|
you should expect). |
1093 |
|
|
|
1094 |
|
|
See C<filter_sir> for more info. |
1095 |
|
|
EOD |
1096 |
|
|
Code => q^ |
1097 |
|
|
int i, n, I; |
1098 |
|
|
$GENERIC() v; |
1099 |
|
|
|
1100 |
|
|
threadloop %{ |
1101 |
|
|
for (n = 0; n < $SIZE(n); n++) |
1102 |
|
|
{ |
1103 |
|
|
v = 0; |
1104 |
|
|
|
1105 |
|
|
loop(an) %{ |
1106 |
|
|
I = n-$a_x(); |
1107 |
|
|
if (I >= 0) |
1108 |
|
|
v += $a_y() * $x(n=>I); |
1109 |
|
|
%} |
1110 |
|
|
|
1111 |
|
|
loop(bn) %{ |
1112 |
|
|
I = n-$b_x(); |
1113 |
|
|
if (I >= 0) |
1114 |
|
|
v -= $b_y() * $y(n=>I); |
1115 |
|
|
%} |
1116 |
|
|
|
1117 |
|
|
$y(n=>n) = v; |
1118 |
|
|
} |
1119 |
|
|
%} |
1120 |
|
|
^ |
1121 |
|
|
; |
1122 |
|
|
|
1123 |
|
|
pp_addpm <<'EOP'; |
1124 |
|
|
|
1125 |
|
|
=head2 filter_fir input, xcoeffs |
1126 |
|
|
|
1127 |
|
|
Apply a fir (finite impulse response) filter to C<input>. This is the same as |
1128 |
|
|
calling: |
1129 |
|
|
|
1130 |
|
|
filter_sir input, xcoeffs, pdl() |
1131 |
|
|
|
1132 |
|
|
=head2 filter_iir input, ycoeffs |
1133 |
|
|
|
1134 |
|
|
Apply a iir (infinite impulse response) filter to C<input>. This is just |
1135 |
|
|
another way of saying: |
1136 |
|
|
|
1137 |
|
|
filter_sir input, pdl(1), ycoeffs |
1138 |
|
|
|
1139 |
|
|
That is, the first member of C<ycoeffs> is being ignored AND SHOULD BE |
1140 |
|
|
SPECIFIED AS ONE FOR FUTURE COMPATIBILITY! |
1141 |
|
|
|
1142 |
|
|
=cut |
1143 |
|
|
|
1144 |
|
|
sub filter_fir($$) { |
1145 |
|
|
filter_sir ($_[0], $_[1], PDL->pdl()); |
1146 |
|
|
} |
1147 |
|
|
|
1148 |
|
|
sub filter_iir($$) { |
1149 |
|
|
filter_sir ($_[0], PDL->pdl(1), $_[1]); |
1150 |
|
|
} |
1151 |
|
|
|
1152 |
|
|
=head2 filter_comb input, delay*, scaler |
1153 |
|
|
|
1154 |
|
|
Apply a comb filter to the piddle C<input>. This is implemented using a |
1155 |
|
|
delay line of length C<delay> (which must be 1 or larger and can be |
1156 |
|
|
non-integer) and a feedback scaler. |
1157 |
|
|
|
1158 |
|
|
y(n) = x(n-size-1) + scaler * y(n-size) |
1159 |
|
|
|
1160 |
|
|
cf. C<filter_notch> and http://www.harmony-central.com/Effects/Articles/Reverb/comb.html |
1161 |
|
|
|
1162 |
|
|
=head2 filter_notch input, delay*, scaler |
1163 |
|
|
|
1164 |
|
|
Apply a comb filter to the piddle C<input>. This is implemented using a |
1165 |
|
|
delay line of length C<delay> (which must be 1 or larger and can be |
1166 |
|
|
non-integer) and a feedforward scaler. |
1167 |
|
|
|
1168 |
|
|
y(n) = x(n-size-1) * scaler + y(n-size) |
1169 |
|
|
|
1170 |
|
|
As a rule of thumb, the decay time of the feedback part is |
1171 |
|
|
C<7*delay/(1-scaler)> samples, so to get a decay of Dur seconds, |
1172 |
|
|
C<scaler <= 1-7*delay/(Dur*Srate)>. The peak gain is C<1/(1-(abs |
1173 |
|
|
scaler))>. The peaks (or valleys in notch's case) are evenly spaced at |
1174 |
|
|
C<srate/delay>. The height (or depth) thereof is determined by scaler |
1175 |
|
|
-- the closer to 1.0, the more pronounced. See Julius Smith's "An |
1176 |
|
|
Introduction to Digital Filter Theory" in Strawn "Digital Audio Signal |
1177 |
|
|
Processing", or Smith's "Music Applications of Digital Waveguides" |
1178 |
|
|
|
1179 |
|
|
=head2 filter_allpass input, delay*, scaler-feedback, scaler-feedforward |
1180 |
|
|
|
1181 |
|
|
C<filter_allpass> or "moving average comb" is just like C<filter_comb> |
1182 |
|
|
but with an added feedforward term. If C<scaler-feedback == 0>, we get a |
1183 |
|
|
moving average comb filter. If both scaler terms == 0, we get a pure delay |
1184 |
|
|
line. |
1185 |
|
|
|
1186 |
|
|
y(n) = feedforward*x(n-1) + x(n-size-1) + feedback*y(n-size) |
1187 |
|
|
|
1188 |
|
|
cf. http://www.harmony-central.com/Effects/Articles/Reverb/allpass.html |
1189 |
|
|
|
1190 |
|
|
=cut |
1191 |
|
|
|
1192 |
|
|
sub _filter_combnotchall { |
1193 |
|
|
my ($pdl, $size, $scalercomb, $scalernotch, $scalerall) = @_; |
1194 |
|
|
my $D = int $size; |
1195 |
|
|
my $F = $size - $D; |
1196 |
|
|
my $a = zeroes $D+2; $a->set($D,1-$F); $a->set($D+1, $F); |
1197 |
|
|
my $b = zeroes $D+2; $b->set($D,$F-1); $b->set($D+1,-$F); |
1198 |
|
|
$a->set(1, $a->at(1) + $scalerall); |
1199 |
|
|
$b->set(0, 1); |
1200 |
|
|
$pdl->filter_sir($a*$scalernotch,$b*$scalercomb); |
1201 |
|
|
} |
1202 |
|
|
|
1203 |
|
|
sub filter_comb { |
1204 |
|
|
_filter_combnotchall $_[0], $_[1], $_[2], 1, 0; |
1205 |
|
|
} |
1206 |
|
|
|
1207 |
|
|
sub filter_notch { |
1208 |
|
|
_filter_combnotchall $_[0], $_[1], 1, $_[2], 0; |
1209 |
|
|
} |
1210 |
|
|
|
1211 |
|
|
sub filter_allpass { |
1212 |
|
|
_filter_combnotchall $_[0], $_[1], $_[2], 1, $_[3]; |
1213 |
|
|
} |
1214 |
|
|
|
1215 |
|
|
=head2 design_remez_fir filter_size, bands(2,b), desired_gain(b), type, [weight(b)] |
1216 |
|
|
|
1217 |
|
|
Calculates the optimal (in the Chebyshev/minimax sense) FIR filter |
1218 |
|
|
impulse response given a set of band edges, the desired reponse on those |
1219 |
|
|
bands, and the weight given to the error in those bands, using the |
1220 |
|
|
Parks-McClellan exchange algorithm. |
1221 |
|
|
|
1222 |
root |
1.6 |
The first argument sets the filter size: C<design_remez_fir> returns as |
1223 |
|
|
many coefficients as specified by this parameter. |
1224 |
root |
1.1 |
|
1225 |
root |
1.6 |
C<bands> is a vector of band edge pairs (start - end), which specify the |
1226 |
root |
1.1 |
start and end of the bands in the filter specification. These must be |
1227 |
|
|
non-overlapping and sorted in increasing order. Only values between C<0> |
1228 |
|
|
(0 Hz) and C<0.5> (the Nyquist frequency) are allowed. |
1229 |
|
|
|
1230 |
|
|
C<des> specifies the desired gain in these bands. |
1231 |
|
|
|
1232 |
|
|
C<weight> can be used to give each band a different weight. If absent, a |
1233 |
|
|
vector of ones is used. |
1234 |
|
|
|
1235 |
|
|
C<type> is any of the exported constants C<BANDPASS>, C<DIFFERENTIATOR> or |
1236 |
|
|
C<HILBERT> and can be used to select various design types (use C<BANDPASS> |
1237 |
|
|
until this is documented ;) |
1238 |
|
|
|
1239 |
|
|
=cut |
1240 |
|
|
|
1241 |
|
|
sub design_remez_fir { |
1242 |
|
|
my ($size, $bands, $des, $type, $weight) = @_; |
1243 |
|
|
$weight = ones $des unless defined $weight; |
1244 |
|
|
$size = zeroes $size; |
1245 |
|
|
_design_remez_fir($bands, $des, $weight, $size, $type); |
1246 |
|
|
$size; |
1247 |
|
|
} |
1248 |
|
|
|
1249 |
|
|
=head2 filter_src input, srate, [width], [sr-mod] |
1250 |
|
|
|
1251 |
|
|
Generic sampling rate conversion, implemented by convoluting C<input> with |
1252 |
|
|
a sinc function of size C<width> (default when unspecified or zero: 5). |
1253 |
|
|
|
1254 |
|
|
C<srate> determines the input rate / output rate ratio, i.e. values > 1 |
1255 |
|
|
speed up, values < 1 slow down. Values < 0 are allowed and reverse the |
1256 |
|
|
signal. |
1257 |
|
|
|
1258 |
|
|
If C<sr_mod> is omitted, the size of the output piddle is calculcated |
1259 |
|
|
as C<length(input)/abs(srate)>, e.g. it provides the full stretched or |
1260 |
|
|
shrinked input signal. |
1261 |
|
|
|
1262 |
|
|
If C<sr_mod> is specified it must be as large as the desired output, i.e. |
1263 |
|
|
it's size determines the output size. Each value in C<sr_mod> is added to |
1264 |
|
|
C<srate> at the given point in "time", so it can be used to "modulate" the |
1265 |
|
|
sampling rate change. |
1266 |
|
|
|
1267 |
|
|
# create a sound effect in the style of "Forbidden Planet" |
1268 |
|
|
$osc = 0.3 * gen_oscil $osc, 30 / $pdl->rate; |
1269 |
|
|
$output = filter_src($input, 1, 0, $osc); |
1270 |
|
|
|
1271 |
|
|
=cut |
1272 |
|
|
|
1273 |
|
|
sub filter_src($$;$$) { |
1274 |
|
|
my ($input, $srate, $width, $sr_mod) = @_; |
1275 |
|
|
|
1276 |
|
|
$width ||= 5; |
1277 |
|
|
$width*2 < $input->getdim(0) or barf "src width too large (> 0.5 * input size)"; |
1278 |
|
|
|
1279 |
|
|
my $output; |
1280 |
|
|
|
1281 |
|
|
if (defined $sr_mod) { |
1282 |
|
|
$output = zeroes $sr_mod; |
1283 |
|
|
} else { |
1284 |
|
|
$output = zeroes ($input->getdim(0) / abs($srate)); |
1285 |
|
|
$sr_mod = PDL->null; |
1286 |
|
|
} |
1287 |
|
|
|
1288 |
|
|
_filter_src($input, $output, $sr_mod, $srate, $width); |
1289 |
|
|
|
1290 |
|
|
$output->rate($input->rate / abs($srate)) if $input->rate; |
1291 |
|
|
$output; |
1292 |
|
|
} |
1293 |
|
|
|
1294 |
|
|
=head2 filter_contrast_enhance input, enhancement |
1295 |
|
|
|
1296 |
|
|
Contrast-enhancement phase-modulates a sound file. It's like audio |
1297 |
|
|
MSG. The actual algorithm is (applied to the normalised sound) |
1298 |
|
|
C<sin(input*pi/2 + (enhancement*sin(input*2*pi)))>. The result is to |
1299 |
|
|
brighten the sound, helping it cut through a huge mix. |
1300 |
|
|
|
1301 |
|
|
=cut |
1302 |
|
|
|
1303 |
|
|
sub filter_contrast_enhance { |
1304 |
|
|
my ($pdl, $idx) = @_; |
1305 |
|
|
$idx = 0.1 unless defined $idx; |
1306 |
|
|
$pdl = $pdl - min $pdl; |
1307 |
|
|
$pdl = ($pdl * (M_PI / max $pdl)) - M_PI * 0.5; |
1308 |
|
|
sin $pdl + $idx * sin $pdl*4; |
1309 |
|
|
} |
1310 |
|
|
|
1311 |
|
|
=head2 filter_granulate input, expansion, [option-hash], option => value... |
1312 |
|
|
|
1313 |
|
|
C<filter_granulate> "granulates" the sound file file. It is the poor |
1314 |
|
|
man's way to change the speed at which things happen in a recorded sound |
1315 |
|
|
without changing the pitches. It works by slicing the input file into |
1316 |
|
|
short pieces, then overlapping these slices to lengthen (or shorten) the |
1317 |
|
|
result; this process is sometimes known as granular synthesis, and is |
1318 |
|
|
similar to the "freeze" function. The duration of each slice is C<length> -- |
1319 |
|
|
the longer, the more like reverb the effect. The portion of the length (on |
1320 |
|
|
a scale from 0 to 1.0) spent on each ramp (up or down) is C<ramp>. This can |
1321 |
|
|
control the smoothness of the result of the overlaps. The more-or-less |
1322 |
|
|
average time between successive segments is C<hop>. The accuracy at which we |
1323 |
|
|
handle this hopping is set by the float C<jitter> -- if C<jitter> is very small, |
1324 |
|
|
you may get an annoying tremolo. The overall amplitude scaler on each |
1325 |
|
|
segment is C<scaler> -- this is used to try to to avoid overflows as we add |
1326 |
|
|
all these zillions of segments together. C<expansion> determines the input |
1327 |
|
|
hop in relation to the output hop; an expansion-amount of C<2.0> should more |
1328 |
|
|
or less double the length of the original, whereas an expansion-amount |
1329 |
|
|
of C<1.0> should return something close to the original speed. |
1330 |
|
|
|
1331 |
|
|
The defaults for the arguments/options are: |
1332 |
|
|
|
1333 |
|
|
expansion 1.0 |
1334 |
|
|
length(*) 0.15 |
1335 |
|
|
scaler 0.6 |
1336 |
|
|
hop(*) 0.05 |
1337 |
|
|
ramp 0.4 |
1338 |
|
|
jitter(*) 0.5 |
1339 |
|
|
maxsize infinity |
1340 |
|
|
|
1341 |
|
|
The parameters/options marked with (*) actually depend on the sampling |
1342 |
|
|
rate, and are always multiplied by the C<rate> attribute of the |
1343 |
|
|
piddle internally. If the piddle lacks that attribute, 44100 is |
1344 |
|
|
assumed. NOTE: This is different to most other filters, but should be ok |
1345 |
|
|
since C<filter_granulate> only makes sense for audiofiles. |
1346 |
|
|
|
1347 |
|
|
=cut |
1348 |
|
|
|
1349 |
|
|
sub filter_granulate { |
1350 |
|
|
my $input = shift; |
1351 |
|
|
my $expansion = shift || 1; |
1352 |
|
|
my %hdr = (); %hdr = (%hdr, %{+shift}) if ref $_[0]; %hdr = (%hdr, @_) if @_; |
1353 |
|
|
|
1354 |
|
|
my $rate = $hdr{rate} || 44100; |
1355 |
|
|
my $length = ($hdr{length} || 0.15) * $rate; |
1356 |
|
|
my $scaler = $hdr{scaler} || 0.6; |
1357 |
|
|
my $hop = ($hdr{hop} || 0.05) * $rate; |
1358 |
|
|
my $ramp = $hdr{ramp} || 0.4; |
1359 |
|
|
my $jitter = ($hdr{jitter} || 0.5) * $rate; |
1360 |
|
|
my $maxsize = $hdr{maxsize}|| 0; |
1361 |
|
|
|
1362 |
|
|
$output = zeroes $expansion * $input->getdim(0); |
1363 |
|
|
|
1364 |
|
|
_filter_granulate($input, $output, |
1365 |
|
|
$expansion, $length, $scaler, |
1366 |
|
|
$hop, $ramp, $jitter, $maxsize); |
1367 |
|
|
$output; |
1368 |
|
|
} |
1369 |
|
|
|
1370 |
|
|
EOP |
1371 |
|
|
|
1372 |
|
|
pp_def '_design_remez_fir', |
1373 |
|
|
Pars => 'bands(2,b); des(b); weight(b); [o]h(n)', |
1374 |
|
|
OtherPars => 'int type', |
1375 |
|
|
GenericTypes => [D], |
1376 |
|
|
Doc => undef, |
1377 |
|
|
PMCode => undef, |
1378 |
|
|
Code => q^ |
1379 |
|
|
int type = $COMP(type); |
1380 |
|
|
|
1381 |
|
|
if (type != BANDPASS |
1382 |
|
|
&& type != DIFFERENTIATOR |
1383 |
|
|
&& type != HILBERT) |
1384 |
|
|
barf ("design_remez_fir: illegal type specified (none of BANDPASS, DIFFERENTIATOR, HILBERT)"); |
1385 |
|
|
|
1386 |
|
|
remez ($P(h), $SIZE(n), $SIZE(b), $P(bands), $P(des), $P(weight), type); |
1387 |
|
|
^ |
1388 |
|
|
; |
1389 |
|
|
|
1390 |
|
|
pp_def '_filter_src', |
1391 |
|
|
Pars => 'input(n); [o]output(m); sr_mod(m2)', |
1392 |
|
|
OtherPars => 'double srate; int width', |
1393 |
|
|
GenericTypes => [D], |
1394 |
|
|
Doc => undef, |
1395 |
|
|
PMCode => undef, |
1396 |
|
|
Code => q! |
1397 |
|
|
mus_src ($P(input), $SIZE(n), $P(output), $SIZE(m), |
1398 |
|
|
$COMP(srate), |
1399 |
|
|
$SIZE(m) == $SIZE(m2) ? $P(sr_mod) : 0, |
1400 |
|
|
$COMP(width)); |
1401 |
|
|
! |
1402 |
|
|
; |
1403 |
|
|
|
1404 |
|
|
pp_def '_filter_granulate', |
1405 |
|
|
Pars => 'input(n); [o]output(m)', |
1406 |
|
|
OtherPars => 'double expansion; double length; double scaler; double hop; double ramp; double jitter; int max_size', |
1407 |
|
|
GenericTypes => [D], |
1408 |
|
|
Doc => undef, |
1409 |
|
|
PMCode => undef, |
1410 |
|
|
Code => q! |
1411 |
|
|
mus_granulate ($P(input), $SIZE(n), $P(output), $SIZE(m), |
1412 |
|
|
$COMP(expansion), $COMP(length), $COMP(scaler), |
1413 |
|
|
$COMP(hop), $COMP(ramp), $COMP(jitter), |
1414 |
|
|
$COMP(max_size)); |
1415 |
|
|
! |
1416 |
|
|
; |
1417 |
|
|
|
1418 |
|
|
pp_addpm <<'EOP'; |
1419 |
|
|
|
1420 |
|
|
=head2 audiomix pos1, data1, pos2, data2, ... |
1421 |
|
|
|
1422 |
|
|
Generate a mix of all given piddles. The resulting piddle will contain the |
1423 |
|
|
sum of all data-piddles at their respective positions, so some scaling |
1424 |
|
|
will be necessary before or after the mixing operation (e.g. scale2short). |
1425 |
|
|
|
1426 |
|
|
# mix the sound gong1 at position 0, the sound bass5 at position 22100 |
1427 |
|
|
# and gong2 at position 44100. The resulting piddle will be large enough |
1428 |
|
|
# to accomodate all the sounds: |
1429 |
|
|
$mix = audiomix 0, $gong1, 44100, $gong2, 22100, $gong2 |
1430 |
|
|
|
1431 |
|
|
=cut |
1432 |
|
|
|
1433 |
|
|
sub audiomix { |
1434 |
|
|
my(@x,@p); |
1435 |
|
|
my($dur,$i,$end,$s,$p); |
1436 |
|
|
for ($i = 0; $i < @_; $i += 2) { |
1437 |
|
|
$end = $_[$i] + $_[$i+1]->getdim(0); |
1438 |
|
|
$dur = $end if $end > $dur; |
1439 |
|
|
} |
1440 |
|
|
my $pdl = zeroes $dur; |
1441 |
|
|
for ($i = 0; $i < @_; $i += 2) { |
1442 |
|
|
($s,$p) = (int($_[$i]), $_[$i+1]); |
1443 |
|
|
($end = $pdl->slice("$s:".($s+$p->getdim(0)-1))) += $p; |
1444 |
|
|
} |
1445 |
|
|
$pdl; |
1446 |
|
|
} |
1447 |
|
|
|
1448 |
|
|
=head2 filter_center piddle |
1449 |
|
|
|
1450 |
|
|
Normalize the piddle so that it is centered around C<y = 0> and has |
1451 |
|
|
maximal amplitude of 1. |
1452 |
|
|
|
1453 |
|
|
=cut |
1454 |
|
|
|
1455 |
|
|
sub filter_center($) { |
1456 |
|
|
my $piddle = shift; |
1457 |
|
|
$piddle = $piddle - min $piddle; |
1458 |
|
|
$piddle * (2 / max $piddle) - 1; |
1459 |
|
|
} |
1460 |
|
|
|
1461 |
|
|
=head2 scale2short piddle |
1462 |
|
|
|
1463 |
|
|
This method takes a sound in any format (preferably float or double) and |
1464 |
|
|
scales it to fit into a signed short value, suitable for playback using |
1465 |
|
|
C<playudio> or similar functions. |
1466 |
|
|
|
1467 |
|
|
=cut |
1468 |
|
|
|
1469 |
|
|
sub scale2short { |
1470 |
|
|
my $pdl = shift->float; |
1471 |
|
|
($pdl * (1 / max abs $pdl) * 32767.5 - 0.5)->short; |
1472 |
|
|
} |
1473 |
|
|
|
1474 |
|
|
=head2 gen_fft_window size*, type, [$beta] |
1475 |
|
|
|
1476 |
|
|
Creates and returns a specific fft window. The C<type> is any of the following. |
1477 |
|
|
These are (case-insensitive) strings, so you might need to quote them. |
1478 |
|
|
|
1479 |
|
|
RECTANGULAR just ones (the identity window) |
1480 |
|
|
HANNING 0.50 - 0.50 * cos (0 .. 2pi) |
1481 |
|
|
HAMMING 0.54 - 0.46 * cos (0 .. 2pi) |
1482 |
|
|
WELCH 1 - (-1 .. 1) ** 2 |
1483 |
|
|
PARZEN the triangle window |
1484 |
|
|
BARTLETT the symmetric triangle window |
1485 |
|
|
BLACKMAN2 blackman-harris window of order 2 |
1486 |
|
|
BLACKMAN3 blackman-harris window of order 3 |
1487 |
|
|
BLACKMAN4 blackman-harris window of order 4 |
1488 |
|
|
EXPONENTIAL the exponential window |
1489 |
|
|
KAISER the kaiser/bessel window (using the parameter C<beta>) |
1490 |
|
|
CAUCHY the cauchy window (using the parameter <beta>) |
1491 |
|
|
POISSON the poisson window (exponential using parameter C<beta>) |
1492 |
|
|
RIEMANN the riemann window (sinc) |
1493 |
|
|
GAUSSIAN the gaussian window of order C<beta>) |
1494 |
|
|
TUKEY the tukey window (C<beta> specifies how much of the window |
1495 |
|
|
consists of ones). |
1496 |
|
|
COSCOS the cosine-squared window (a partition of unity) |
1497 |
|
|
SINC same as RIEMANN |
1498 |
|
|
HANN same as HANNING (his name was Hann, not Hanning) |
1499 |
|
|
|
1500 |
|
|
LIST this "type" is special in that it returns a list of all types |
1501 |
|
|
|
1502 |
|
|
=cut |
1503 |
|
|
|
1504 |
|
|
sub gen_fft_window { |
1505 |
|
|
my ($size, $type, $beta) = @_; |
1506 |
|
|
|
1507 |
|
|
$beta = 2.5 unless defined $beta; |
1508 |
|
|
|
1509 |
|
|
$size = $size->getdim(0) if ref $size; |
1510 |
|
|
$size > 2 or barf "fft window size too small"; |
1511 |
|
|
|
1512 |
root |
1.5 |
my $midn = $size >> 1; |
1513 |
root |
1.1 |
my $midm1 = ($size-1) >> 1; |
1514 |
|
|
my $midp1 = ($size+1) >> 1; |
1515 |
root |
1.5 |
my $dur = zeroes $size; |
1516 |
root |
1.1 |
my $sf = ($size-1)/$size; |
1517 |
|
|
%fft_window = ( |
1518 |
|
|
RECTANGULAR => sub { |
1519 |
|
|
$dur->ones |
1520 |
|
|
}, |
1521 |
|
|
HANNING => sub { |
1522 |
|
|
0.5 - 0.5 * cos $dur->xlinvals(0,M_2PI*$sf) |
1523 |
|
|
}, |
1524 |
|
|
HAMMING => sub { |
1525 |
|
|
0.53836 - 0.46164 * cos $dur->xlinvals(0,M_2PI*$sf) |
1526 |
|
|
}, |
1527 |
|
|
WELCH => sub { |
1528 |
|
|
1 - $dur->xlinvals(-1,$sf)**2; |
1529 |
|
|
}, |
1530 |
|
|
PARZEN => sub { |
1531 |
|
|
1 - abs $dur->xlinvals(-$sf,$sf); |
1532 |
|
|
}, |
1533 |
|
|
BARTLETT => sub { |
1534 |
|
|
1 - abs $dur->xlinvals(-1,1); |
1535 |
|
|
}, |
1536 |
|
|
BLACKMAN2 => sub { |
1537 |
|
|
my $cx = cos $dur->xlinvals(0,M_2PI*$sf); |
1538 |
|
|
0.34401 + ($cx * (-0.49755 + ($cx * 0.15844))); |
1539 |
|
|
}, |
1540 |
|
|
BLACKMAN3 => sub { |
1541 |
|
|
my $cx = cos $dur->xlinvals(0,M_2PI*$sf); |
1542 |
|
|
0.21747 + ($cx * (-0.45325 + ($cx * (0.28256 - ($cx * 0.04672))))); |
1543 |
|
|
}, |
1544 |
|
|
BLACKMAN4 => sub { |
1545 |
|
|
my $cx = cos $dur->xlinvals(0,M_2PI*$sf); |
1546 |
|
|
0.084037 + ($cx * (-0.29145 + ($cx * (0.375696 + ($cx * (-0.20762 + ($cx * 0.041194))))))); |
1547 |
|
|
}, |
1548 |
|
|
EXPONENTIAL => sub { |
1549 |
|
|
2 ** (1 - abs $dur->xlinvals(-1,$sf)) - 1; |
1550 |
|
|
}, |
1551 |
|
|
KAISER => sub { |
1552 |
|
|
bessi0 ($beta * sqrt(1 - $dur->xlinvals(-1,$sf)**2)) / bessi0 ($beta); |
1553 |
|
|
}, |
1554 |
|
|
CAUCHY => sub { |
1555 |
|
|
1 / (1 + ($dur->xlinvals(-1,$sf)*$beta)**2); |
1556 |
|
|
}, |
1557 |
|
|
POISSON => sub { |
1558 |
|
|
exp (-$beta * abs $dur->xlinvals(-1,$sf)); |
1559 |
|
|
}, |
1560 |
|
|
RIEMANN => sub { |
1561 |
|
|
my $dur1 = $dur->slice("0:$midm1"); |
1562 |
|
|
my $dur2 = $dur->slice("-1:$midp1:-1"); |
1563 |
|
|
my $cx; |
1564 |
|
|
$cx = $dur1->xlinvals(M_PI,M_2PI/$size); $dur1 .= sin ($cx) / $cx; |
1565 |
|
|
$cx = $dur2->xlinvals(M_PI,M_2PI/$size); $dur2 .= sin ($cx) / $cx; |
1566 |
|
|
($cx = $dur->slice("$midm1:$midp1")) .= 1; |
1567 |
|
|
$dur; |
1568 |
|
|
}, |
1569 |
|
|
GAUSSIAN => sub { |
1570 |
|
|
exp (-0.5 * ($beta * abs $dur->xlinvals(-1,$sf))**2); |
1571 |
|
|
}, |
1572 |
|
|
TUKEY => sub { |
1573 |
|
|
$beta >= 0 && $beta <= 1 or barf "beta must be between 0 and 1 for the tukey window"; |
1574 |
|
|
my $cx = int ($midn * (1 - $beta)); |
1575 |
|
|
my $cX = $size-$cx-1; |
1576 |
|
|
my $dur1 = $dur->slice("0:".($cx-1)); |
1577 |
|
|
my $dur2 = $dur->slice("-1:".($cX+1).":-1"); |
1578 |
|
|
my $dur3 = $dur->slice("$cx:$cX"); |
1579 |
|
|
|
1580 |
|
|
$dur1 .= 0.5 * (1 - cos $dur1->xlinvals(0, M_PI)); |
1581 |
|
|
$dur2 .= 0.5 * (1 - cos $dur2->xlinvals(0, M_PI)); |
1582 |
|
|
$dur3 .= 1; |
1583 |
|
|
|
1584 |
|
|
$dur; |
1585 |
|
|
}, |
1586 |
|
|
COSCOS => sub { |
1587 |
|
|
(cos $dur->xlinvals(M_PI/-2,M_PI/2))**2; |
1588 |
|
|
}, |
1589 |
|
|
LIST => sub { |
1590 |
|
|
grep $_ ne "LIST", keys %fft_window; |
1591 |
|
|
}, |
1592 |
|
|
); |
1593 |
|
|
|
1594 |
|
|
$fft_window{SINC} = $fft_window{RIEMANN}; |
1595 |
|
|
$fft_window{HANN} = $fft_window{HANNING}; |
1596 |
|
|
|
1597 |
|
|
$type = uc $type; |
1598 |
|
|
$fft_window{$type} or barf "$type: no such window"; |
1599 |
|
|
$fft_window{$type}->(); |
1600 |
|
|
} |
1601 |
|
|
|
1602 |
|
|
=head2 cplx(2,n) = rfft real(n) |
1603 |
|
|
|
1604 |
|
|
Do a (complex fft) of C<real> (extended to complex so that the imaginary |
1605 |
|
|
part is zero), and return the complex fft result. This function tries |
1606 |
|
|
to use L<PDL::FFTW|PDL::FFTW> (which is faster for large vectors) when |
1607 |
|
|
available, and falls back to L<PDL::FFT|PDL::FFT>, which is likely |
1608 |
|
|
to return different phase signs (due to different kernel functions), |
1609 |
|
|
so beware! In fact, since C<rfft> has to shuffle the data when using |
1610 |
|
|
PDL::FFT, the fallback is always slower. |
1611 |
|
|
|
1612 |
|
|
When using PDL::FFTW, a wisdom file ~/.pdl_wisdom is used and updated, if |
1613 |
|
|
possible. |
1614 |
|
|
|
1615 |
|
|
=head2 real(n) = irfft cplx(2,n) |
1616 |
|
|
|
1617 |
|
|
The inverse transformation (see C<rfft>). C<irfft rfft $pdl == $pdl> |
1618 |
|
|
always holds. |
1619 |
|
|
|
1620 |
|
|
=cut |
1621 |
|
|
|
1622 |
|
|
my $fftw_loaded; |
1623 |
|
|
sub _fftw { |
1624 |
|
|
defined $fftw_loaded or eval { |
1625 |
|
|
$fftw_loaded = 0; |
1626 |
|
|
require PDL::FFTW; |
1627 |
root |
1.2 |
PDL::FFTW::load_wisdom("$ENV{HOME}/.pdl_wisdom") |
1628 |
|
|
if -r "$ENV{HOME}/.pdl_wisdom"; |
1629 |
root |
1.1 |
$fftw_loaded = 1; |
1630 |
|
|
}; |
1631 |
|
|
$fftw_loaded; |
1632 |
|
|
} |
1633 |
|
|
|
1634 |
|
|
sub rfft { |
1635 |
|
|
my $data = $_[0]; |
1636 |
|
|
if (_fftw) { |
1637 |
|
|
my $x = $data->r2C; |
1638 |
root |
1.3 |
$x = PDL::FFTW::fftw $x; |
1639 |
root |
1.1 |
$x; |
1640 |
|
|
} else { |
1641 |
|
|
require PDL::FFT; |
1642 |
|
|
my @fft = ($data->copy, $data->zeroes); |
1643 |
|
|
PDL::FFT::fft(@fft); |
1644 |
|
|
cat(@fft)->xchg(0,1); |
1645 |
|
|
} |
1646 |
|
|
} |
1647 |
|
|
|
1648 |
|
|
sub irfft { |
1649 |
|
|
if (_fftw) { |
1650 |
|
|
$x = $_[0]->copy; |
1651 |
root |
1.3 |
$x = PDL::FFTW::ifftw $x; |
1652 |
root |
1.1 |
re $x / $x->getdim(1); |
1653 |
|
|
} else { |
1654 |
|
|
require PDL::FFT; |
1655 |
|
|
my @fft = $_[0]->xchg(0,1)->dog({Break => 1}); |
1656 |
|
|
PDL::FFT::ifft(@fft); |
1657 |
|
|
$fft[0]; |
1658 |
|
|
} |
1659 |
|
|
} |
1660 |
|
|
|
1661 |
|
|
=head2 spectrum data, [norm], [window], [beta] |
1662 |
|
|
|
1663 |
|
|
Returns the spectrum of a given pdl. If C<norm> is absent (or C<undef>), |
1664 |
|
|
it returns the magnitude of the fft of C<data>. When C<norm> == 1 (or |
1665 |
root |
1.6 |
C<eq 'NORM'>, case-insensitive), it returns the magnitude, normalized to be |
1666 |
|
|
between zero and one. If C<norm> == 0 (or C<eq 'dB'>, case-insensitive), then |
1667 |
root |
1.1 |
it returns the magnitude in dB. |
1668 |
|
|
|
1669 |
|
|
C<data> is multiplied with C<window> (if not C<undef>) before calculating |
1670 |
|
|
the fft, and usually contains a window created with C<gen_fft_window> |
1671 |
|
|
(using C<beta>). If C<window> is a string, it is handed over to |
1672 |
|
|
C<gen_fft_window> (together with the beta parameter) to create a window of |
1673 |
|
|
suitable size. |
1674 |
|
|
|
1675 |
|
|
This function could be slightly faster. |
1676 |
|
|
|
1677 |
|
|
=cut |
1678 |
|
|
|
1679 |
|
|
sub spectrum { |
1680 |
|
|
my ($data, $norm, $window, $beta) = @_; |
1681 |
|
|
my $len; |
1682 |
|
|
if (defined $window) { |
1683 |
root |
1.4 |
$window = gen_fft_window ($data->getdim (0), $window, $beta) unless ref $window; |
1684 |
root |
1.1 |
$data = $data * $window; |
1685 |
root |
1.4 |
$len = $window->getdim (0); |
1686 |
root |
1.1 |
} else { |
1687 |
root |
1.4 |
$len = $data->getdim (0); |
1688 |
root |
1.1 |
} |
1689 |
root |
1.4 |
$data = rfft ( |
1690 |
|
|
$data->slice ("0:" . ($len - 1)) |
1691 |
|
|
->sever |
1692 |
|
|
) |
1693 |
|
|
->slice (",0:" . int ($len / 2)) |
1694 |
root |
1.5 |
->PDL::Complex::Cr2p |
1695 |
|
|
->slice ("(0)"); |
1696 |
root |
1.1 |
if ($norm == 1 || lc $norm eq "norm") { |
1697 |
|
|
$data / max $data; |
1698 |
|
|
} elsif (($norm =~ /^[.0]+$/) || (lc $norm eq "db")) { |
1699 |
|
|
log (1e-37 + $data / max $data) * (20 / log 10); |
1700 |
|
|
} else { |
1701 |
|
|
$data; |
1702 |
|
|
} |
1703 |
|
|
} |
1704 |
|
|
|
1705 |
|
|
=head2 concat pdl, pdl... |
1706 |
|
|
|
1707 |
|
|
This is not really an audio-related function. It simply takes all piddles |
1708 |
|
|
and concats them into a larger one. At the moment it only supports |
1709 |
|
|
single-dimensional piddles and is implemented quite slowly using perl and |
1710 |
|
|
data-copying, but that might change... |
1711 |
|
|
|
1712 |
|
|
=cut |
1713 |
|
|
|
1714 |
|
|
sub concat { |
1715 |
|
|
my $len = sum pdl map $_->getdim(0), @_; |
1716 |
|
|
my $r = zeroes $len; |
1717 |
|
|
my $o = 0; |
1718 |
|
|
for (@_) { |
1719 |
|
|
my $e = $o + $_->getdim(0) - 1; |
1720 |
|
|
(my $t = $r->slice("$o:$e")) .= $_; |
1721 |
|
|
$o = $e + 1; |
1722 |
|
|
} |
1723 |
|
|
$r; |
1724 |
|
|
} |
1725 |
|
|
|
1726 |
|
|
EOP |
1727 |
|
|
|
1728 |
|
|
# TODO: document! |
1729 |
|
|
pp_def 'filter_convolve', |
1730 |
|
|
Pars => 'input(n); kernel(m); int fftsize(); [o]output(n)', |
1731 |
|
|
GenericTypes => [D], |
1732 |
|
|
Code => q! |
1733 |
|
|
mus_convolve ($P(input), $P(output), $SIZE(n), $P(kernel), $fftsize(), $SIZE(m)); |
1734 |
|
|
! |
1735 |
|
|
; |
1736 |
|
|
|
1737 |
|
|
pp_def 'rshift', |
1738 |
|
|
Doc => <<'EOD', |
1739 |
|
|
=for ref |
1740 |
|
|
|
1741 |
|
|
Shift vector elements without wrap and fill the free space with a |
1742 |
|
|
constant. Flows data back & forth, for values that overlap. |
1743 |
|
|
|
1744 |
|
|
Positive values shift right, negative values shift left. |
1745 |
|
|
|
1746 |
|
|
EOD |
1747 |
|
|
Pars => 'x(n); int shift(); c(); [oca]y(n)', |
1748 |
|
|
DefaultFlow => 1, |
1749 |
|
|
Reversible => 1, |
1750 |
|
|
PMCode => ' |
1751 |
|
|
sub PDL::rshift { |
1752 |
|
|
my @a = @_; |
1753 |
|
|
if($#a == 3) { |
1754 |
|
|
&PDL::_rshift_int(@a);@a=(); |
1755 |
|
|
} elsif($#a == 1 || $#a == 2) { |
1756 |
|
|
$a[3] = PDL->nullcreate($a[0]); |
1757 |
|
|
&PDL::_rshift_int(@a); |
1758 |
|
|
$a[3]; |
1759 |
|
|
} else { |
1760 |
|
|
barf "Invalid number of arguments for shiftin"; |
1761 |
|
|
} |
1762 |
|
|
} |
1763 |
|
|
', |
1764 |
|
|
Code=>' |
1765 |
|
|
int i,j; |
1766 |
|
|
int n_size = $SIZE(n); |
1767 |
|
|
for(i = -$shift(), j=0; j < n_size; i++, j++) |
1768 |
|
|
$y(n=>j) = i >= 0 && i < n_size ? $x(n=>i) : $c(); |
1769 |
|
|
', |
1770 |
|
|
BackCode=>' |
1771 |
|
|
int i,j; |
1772 |
|
|
int n_size = $SIZE(n); |
1773 |
|
|
for(i = -$shift(), j=0; j < n_size; i++, j++) |
1774 |
|
|
if (i >= 0 && i < n_size) |
1775 |
|
|
$x(n=>i) = $y(n=>j); |
1776 |
|
|
' |
1777 |
|
|
; |
1778 |
|
|
|
1779 |
|
|
pp_def 'polynomial', |
1780 |
|
|
Pars => 'coeffs(n); x(m); [o]out(m)', |
1781 |
|
|
Doc => 'evaluate the polynomial with coefficients C<coeffs> at the position(s) C<x>. C<coeffs[0]> is the constant term.', |
1782 |
|
|
Code => q! |
1783 |
|
|
loop(m) %{ |
1784 |
|
|
$GENERIC() x = 1; |
1785 |
|
|
$GENERIC() o = 0; |
1786 |
|
|
|
1787 |
|
|
loop(n) %{ |
1788 |
|
|
o += $coeffs() * x; |
1789 |
|
|
x *= $x(); |
1790 |
|
|
%} |
1791 |
|
|
|
1792 |
|
|
$out() = o; |
1793 |
|
|
%} |
1794 |
|
|
! |
1795 |
|
|
; |
1796 |
|
|
|
1797 |
|
|
pp_def 'linear_interpolate', |
1798 |
|
|
Pars => 'x(); fx(n); fy(n); [o]y()', |
1799 |
|
|
GenericTypes => [L,F,D], |
1800 |
|
|
Doc => ' |
1801 |
|
|
Look up the ordinate C<x> in the function given by C<fx> and C<fy> |
1802 |
|
|
and return a linearly interpolated value (somewhat optimized for many |
1803 |
|
|
lookups). |
1804 |
|
|
|
1805 |
|
|
C<fx> specifies the ordinates (x-coordinates) of the function and most be |
1806 |
|
|
sorted in increasing order. C<fy> are the y-coordinates of the function in |
1807 |
|
|
these points. |
1808 |
|
|
', |
1809 |
|
|
Code => q! |
1810 |
|
|
int d = 0; |
1811 |
|
|
int D = $SIZE(n) - 1; /* a tribute to PP's stupidity */ |
1812 |
|
|
$GENERIC() xmin =$fx(n=>0); |
1813 |
|
|
$GENERIC() xmax =$fx(n=>D); |
1814 |
|
|
|
1815 |
|
|
threadloop %{ |
1816 |
|
|
$GENERIC() x = $x(); |
1817 |
|
|
$GENERIC() x1, x2, y1, y2; |
1818 |
|
|
|
1819 |
|
|
if (x <= xmin) |
1820 |
|
|
$y() = $fy(n=>0); |
1821 |
|
|
else if (x >= xmax) |
1822 |
|
|
$y() = $fy(n=>D); |
1823 |
|
|
else |
1824 |
|
|
{ |
1825 |
|
|
while ($fx(n=>d) > x) d--; |
1826 |
|
|
while ($fx(n=>d) <= x) d++; |
1827 |
|
|
/* 0 < d <= D */ |
1828 |
|
|
/* fx(d-1) < x <= fx(d) */ |
1829 |
|
|
x2 = $fx(n=>d); |
1830 |
|
|
y2 = $fy(n=>d); |
1831 |
|
|
d--; |
1832 |
|
|
x1 = $fx(n=>d); |
1833 |
|
|
y1 = $fy(n=>d); |
1834 |
|
|
$y() = y1 + (x-x1)*(y2-y1)/(x2-x1); |
1835 |
|
|
} |
1836 |
|
|
%} |
1837 |
|
|
! |
1838 |
|
|
; |
1839 |
|
|
|
1840 |
|
|
pp_def 'bessi0', |
1841 |
|
|
Pars => 'a(); [o]b()', |
1842 |
|
|
Doc => 'calculate the (approximate) modified bessel function of the first kind', |
1843 |
|
|
GenericTypes => [F,D], |
1844 |
|
|
Code => q! |
1845 |
|
|
$GENERIC() x = $a(); |
1846 |
|
|
if (x > -3.75 && x < 3.75) |
1847 |
|
|
{ |
1848 |
|
|
$GENERIC() y = (x / 3.75); y *= y; |
1849 |
|
|
$b() = 1 + (y * (3.5156229 + |
1850 |
|
|
(y * (3.0899424 + |
1851 |
|
|
(y * (1.2067492 + |
1852 |
|
|
(y * (0.2659732 + |
1853 |
|
|
(y * (0.360768e-1 + |
1854 |
|
|
(y * 0.45813e-2))))))))))); |
1855 |
|
|
} |
1856 |
|
|
else |
1857 |
|
|
{ |
1858 |
|
|
$GENERIC() ax = x < 0 ? -x : x; |
1859 |
|
|
$GENERIC() y = ax / 3.75; |
1860 |
|
|
$b() = ((exp(ax) / sqrt(ax)) * |
1861 |
|
|
(0.39894228 + |
1862 |
|
|
(y * (0.1328592e-1 + |
1863 |
|
|
(y * (0.225319e-2 + |
1864 |
|
|
(y * (-0.157565e-2 + |
1865 |
|
|
(y * (0.916281e-2 + |
1866 |
|
|
(y * (-0.2057706e-1 + |
1867 |
|
|
(y * (0.2635537e-1 + |
1868 |
|
|
(y * (-0.1647633e-1 + |
1869 |
|
|
(y * 0.392377e-2))))))))))))))))); |
1870 |
|
|
} |
1871 |
|
|
!, |
1872 |
|
|
; |
1873 |
|
|
|
1874 |
|
|
pp_def 'fast_sin', |
1875 |
|
|
Pars => 'r(n); [o]s(n)', |
1876 |
|
|
GenericTypes => [F,D], |
1877 |
|
|
Doc => 'fast sine function (inaccurate table lookup with ~12 bits precision)', |
1878 |
|
|
Code => q^ |
1879 |
|
|
# define SINE_SIZE 16384 |
1880 |
|
|
static float *sine_table = 0; |
1881 |
|
|
|
1882 |
|
|
if (!sine_table) |
1883 |
|
|
{ |
1884 |
|
|
int i; |
1885 |
|
|
Float phase, incr; |
1886 |
|
|
sine_table = (float *) calloc (SINE_SIZE + 1, sizeof (float)); |
1887 |
|
|
incr = M_2PI / (Float) SINE_SIZE; |
1888 |
|
|
for (i = 0, phase = 0.0; i < SINE_SIZE + 1; i++, phase += incr) |
1889 |
|
|
sine_table[i] = (float) sin (phase); |
1890 |
|
|
} |
1891 |
|
|
|
1892 |
|
|
threadloop %{ |
1893 |
|
|
loop(n) %{ |
1894 |
|
|
$s() = sine_table[((int)($r() * (SINE_SIZE / M_2PI)) % SINE_SIZE + SINE_SIZE) % SINE_SIZE]; |
1895 |
|
|
%} |
1896 |
|
|
%}^ |
1897 |
|
|
; |
1898 |
|
|
|
1899 |
|
|
pp_addpm {At => Bot}, <<'EOD'; |
1900 |
|
|
|
1901 |
|
|
=head1 AUTHOR |
1902 |
|
|
|
1903 |
root |
1.7 |
Marc Lehmann <schmorp@schmorp.de>. The ideas were mostly taken from common |
1904 |
root |
1.1 |
lisp music (CLM), by Bill Schottstaedt C<bil@ccrma.stanford.edu>. I also |
1905 |
|
|
borrowed many explanations (and references) from the clm docs and some |
1906 |
|
|
code from clm.c. Highly inspiring! |
1907 |
|
|
|
1908 |
|
|
=head1 SEE ALSO |
1909 |
|
|
|
1910 |
|
|
perl(1), L<PDL>. |
1911 |
|
|
|
1912 |
|
|
=cut |
1913 |
|
|
EOD |
1914 |
|
|
|
1915 |
|
|
pp_addxs "\nINCLUDE: xlib.xs\n"; |
1916 |
|
|
|
1917 |
|
|
pp_done; |
1918 |
|
|
|