1# dsp.rb -- dsp.scm --> dsp.rb 2 3# Translator: Michael Scholz <mi-scholz@users.sourceforge.net> 4# Created: 2005/03/07 13:50:44 5# Changed: 2018/04/15 22:52:00 6 7# comments are taken mostly from dsp.scm 8# 9# module Dsp 10# src_duration(e) 11# dolph(n, gamma) 12# dolph_1(n, gamma) 13# down_oct(n, snd, chn) 14# edot_product(freq, data) 15# stretch_sound_via_dft(factor, snd, chn) 16# compute_uniform_circular_string(size, x0, x1, x2, mass, xspring, damp) 17# testunif(mass, xspring, damp) 18# test_scanned_synthesis(amp, dur, mass, xspring, damp) 19# compute_string(size, x0, x1, x2, masses, xsprings, esprings, damps, haptics) 20# freqdiv(n, snd, chn) 21# adsat(size, beg, dur, snd, chn) 22# spike(snd, chn) 23# spot_freq(samp, snd, chn) 24 25# class Flanger < Musgen 26# initialize(time, amount, speed) 27# inspect 28# to_s 29# run_func(val1, val2) 30# flanger(inval) 31 32# make_flanger(time, amount, speed) 33# flanger?(obj) 34# flanger(gen, inval) 35 36# chorus(size) 37# chordalize(amount, base, chord) 38 39# zero_phase(snd, chn) 40# rotate_phase(func, snd, chn) 41 42# class Asyfm < Musgen 43# initialize(*args) 44# inspect 45# to_s 46# run_func(val1, val2) 47# asyfm_J(input) 48# asyfm_I(input) 49 50# make_asyfm(*args) 51# asyfm?(obj) 52# asyfm_J(gen, input) 53# asyfm_I(gen, input) 54 55# cosine_summation(gen, r) 56# kosine_summation(gen, r, k) 57# fejer_sum(angle, n) 58# legendre_sum(angle, n) 59 60# sum_of_n_sines(angle, n) 61# sum_of_n_odd_sines(angle, n) 62# sum_of_n_odd_cosines(angle, n) 63 64# band_limited_sawtooth(x, a, n, fi) 65# brighten_slightly(amount, snd, chn) 66# brighten_slightly_1(coeffs, snd, chn) 67# spectrum2coeffs(order, spectr) 68# fltit_1(order, spectr) 69 70# make_hilbert_transform(len) 71# make_highpass(fc, len) 72# make_lowpass(fc, len) 73# make_bandpass(flo, fhi, len) 74# make_bandstop(flo, fhi, len) 75# make_differentiator(len) 76 77# make_butter_high_pass(freq) 78# make_butter_low_pass(freq) 79# make_butter_band_pass(freq, band) 80# make_butter_band_reject(freq, band) 81 82# make_biquad(a0, a1, a2, b1, b2) 83# make_iir_low_pass_2(fc, din) 84# make_iir_high_pass_2(fc, din) 85# make_iir_band_pass_2(f1, f2) 86# make_iir_band_stop_2(f1, f2) 87 88# make_eliminate_hum(hum_freq, hum_harmonics, bandwidth) 89# eliminate_hum(gens, x0) 90# make_peaking_2(f1, f2, m) 91# cascade2canonical(a) 92 93# make_butter_lp(m, fc) 94# make_butter_hp(m, fc) 95# make_butter_bp(m, f1, f2) 96# make_butter_bs(m, f1, f2) 97 98# make_notch_frequency_response(cur_srate, freqs, notch_width) 99# notch_channel(freqs, filter_order, beg, dur, snd, chn, edpos, trc, n_width) 100# notch_sound(freqs, filter_order, snd, chn, notch_width) 101# notch_selection(freqs, filter_order, snd, chn, notch_width) 102# fractional_fourier_transform(fr, fi, n, v) 103# z_transform(f, n, z) 104# dht(data) 105# find_sine(freq, beg, dur, snd) 106# goertzel(freq, beg, dur, snd) 107# make_spencer_filter 108 109# any_random(amount, e) 110# gaussian_distribution(s) 111# pareto_distribution(a) 112# inverse_integrate(dist, data_size, e_size) 113# gaussian_envelope(s) 114 115# channel_mean(snd, chn) 116# channel_total_energy(snd, chn) 117# channel_average_power(snd, chn) 118# channel_rms(snd, chn) 119# channel_variance(snd, chn) 120# channel_norm(snd, chn) 121# channel_lp(p, snd, chn) 122# channel_lp_inf(snd, chn) 123# channel2_inner_product(s1, c1, s2, c2) 124# channel2_angle(s1, c1, s2, c2) 125# channel2_orthogonal?(s1, c1, s2, c2) 126# channel2_coefficient_of_projection(s1, c1, s2, c2) 127# channel_distance(s1, c1, s2, c2) 128 129# periodogram(n) 130# shift_channel_pitch(freq, order, beg, dur, snd, chn, edpos) 131# hz_to_2pi(freq) 132# ssb_bank(old_freq, new_freq, pairs, order, bw, beg, dur, snd, chn, edpos) 133# ssb_bank_env(old_freq, new_freq, fq_env, pairs, ord, bw, beg, dur, s, c, ep) 134# 135# vct_polynomial(v, coeffs) 136# channel_polynomial(coeffs, snd, chn) 137# spectral_polynomial(coeffs, snd, chn) 138# scentroid(file, beg, dur, db_floor, rfreq, fftsize) 139# invert_filter(fcoeffs) 140# 141# class Volterra_filter < Musgen 142# initialize(acoeffs, bcoeffs) 143# inspect 144# to_s 145# run_func(val1, val2) 146# volterra_filter(x) 147# 148# make_volterra_filter(acoeffs, bcoeffs) 149# volterra_filter(flt, x) 150# make_moving_sum(size) 151# moving_sum(gen, y) 152# make_moving_rms(size) 153# moving_rms(gen, y) 154# make_moving_length(size) 155# moving_length(gen, y) 156# harmonicizer(freq, coeffs, pairs, order, bw, beg, dur, snd, chn, edpos) 157# linear_src_channel(srinc, snd, chn) 158# 159# class Mfilter < Musgen 160# initialize(decay, freq) 161# inspect 162# to_s 163# mfilter(x_input, y_input) 164# 165# make_mfilter(*args) 166# mfilter(m, x_input, y_input) 167# 168# class Display_bark_fft 169# display_bark_fft(snd, chn) 170# mark_bark_labels(snd, chn) 171# choose_bark_ticks(snd, chn, button, state, x, y, axis) 172# 173# display_bark_fft(off) 174# undisplay_bark_fft 175 176require "ws" 177require "env" 178 179module Dsp 180 # src_duration (see src-channel in extsnd.html) 181 add_help(:src_duration, 182 "src_duration(envelope) \ 183Returns the new duration of a sound after using ENVELOPE for \ 184time-varying sampling-rate conversion.") 185 def src_duration(e) 186 e.map! do |x| 187 x.to_f 188 end 189 ex0 = e.first 190 ex1 = e[-2] 191 all_x = ex1 - ex0 192 dur = 0.0 193 0.step(e.length - 3, 2) do |i| 194 x0, xy0, x1, xy1 = e[i, 4] 195 y0 = xy0.zero? ? 1.0 : (1.0 / xy0) 196 y1 = xy1.zero? ? 1.0 : (1.0 / xy1) 197 area = if (xy0 - xy1).abs < 0.0001 198 y0 * ((x1 - x0) / all_x) 199 else 200 ((log(y1) - log(y0)) / (xy0 - xy1)) * ((x1 - x0) / all_x) 201 end 202 dur += area.abs 203 end 204 dur 205 end 206 207 # Dolph-Chebyshev window 208 # 209 # formula taken from Richard Lyons, "Understanding DSP" 210 # see clm.c for C version (using either GSL's or GCC's complex trig functions) 211 add_help(:dolph, 212 "dolph(n, gamma) \ 213Produces a Dolph-Chebyshev FFT data window of N points \ 214using GAMMA as the window parameter.") 215 def dolph(n, gamma) 216 alpha = cosh(acosh(10.0 ** gamma) / n) 217 den = 1.0 / cosh(n * acosh(alpha)) 218 freq = PI / n 219 rl = make_vct(n) 220 im = make_vct(n) 221 phase = 0.0 222 n.times do |i| 223 val = den * cos(n * acos(alpha * cos(phase))) 224 rl[i] = val.real 225 im[i] = val.imag 226 phase += freq 227 end 228 fft(rl, im, -1) 229 vct_scale!(rl, 1.0 / vct_peak(rl)) 230 j = n / 2 231 n.times do |i| 232 im[i] = rl[j] 233 j += 1 234 if j == n 235 j = 0 236 end 237 end 238 im 239 end if defined? acosh 240 241 # this version taken from Julius Smith's "Spectral Audio..." with 242 # three changes it does the DFT by hand, and is independent of 243 # anything from Snd (fft, vcts etc) 244 add_help(:dolph_1, 245 "dolph_1(n, gamma) \ 246Produces a Dolph-Chebyshev FFT data window of N points \ 247using GAMMA as the window parameter.") 248 def dolph_1(n, gamma) 249 alpha = cosh(acosh(10.0 ** gamma) / n) 250 den = 1.0 / cosh(n * acosh(alpha)) 251 freq = PI / n 252 vals = make_array(n) 253 w = make_array(n) 254 pk = 0.0 255 mult = -1.0 256 phase = -HALF_PI 257 n.times do |i| 258 vals[i] = mult * den * cos(n * acos(alpha * cos(phase))) 259 mult *= -1.0 260 phase += freq 261 end 262 n.times do |i| 263 sum = 0.0 264 n.times do |j| 265 sum = sum + vals[j] * exp((2.0 * Complex(0.0, 1.0) * PI * j * i) / n) 266 end 267 w[i] = sum.abs 268 if w[i] > pk 269 pk = w[i] 270 end 271 end 272 w.map! do |val| 273 val / pk 274 end 275 end if defined? acosh 276 277 # move sound down by n (a power of 2) 278 # I think this is "stretch" in DSP jargon -- to interpolate in the 279 # time domain we're squeezing the frequency domain the power-of-2 280 # limitation is based on the underlying fft function's insistence on 281 # power-of-2 data sizes see stretch-sound-via-dft below for a 282 # general version 283 add_help(:down_oct, 284 "down_oct(n, snd=false, chn=false) \ 285Moves a sound down by power of 2 N.") 286 def down_oct(n, snd = false, chn = false) 287 len = framples(snd, chn) 288 pow2 = (log(len) / log(2)).ceil 289 fftlen = (2 ** pow2).round 290 fftscale = 1.0 / fftlen 291 rl1 = channel2vct(0, fftlen, snd, chn) 292 im1 = make_vct(fftlen) 293 fft(rl1, im1, 1) 294 vct_scale!(rl1, fftscale) 295 vct_scale!(im1, fftscale) 296 rl2 = make_vct(2 * fftlen) 297 im2 = make_vct(2 * fftlen) 298 k = fftlen - 1 299 j = fftlen * n - 1 300 (0...(fftlen / 2)).each do |i| 301 vct_set!(rl2, i, rl1[i]) 302 vct_set!(rl2, j, rl1[k]) 303 vct_set!(im2, i, im1[i]) 304 vct_set!(im2, j, im1[k]) 305 k -= 1 306 j -= 1 307 end 308 fft(rl2, im2, -1) 309 vct2channel(rl2, 0, n * len, snd, chn, false, 310 format("%s(%s", get_func_name, n)) 311 end 312 313 add_help(:edot_product, 314 "edot_product(freq, data) \ 315Sum of (e^freq*i) * data[i]") 316 def edot_product(freq, data) 317 sum = 0.0 318 data.each_with_index do |val, i| 319 sum = sum + exp(i.to_f * freq) * val 320 end 321 sum 322 end unless defined? edot_product 323 324 add_help(:stretch_sound_via_dft, 325 "stretch_sound_via_dft(factor, snd=false, chn=false) \ 326Makes the given channel longer (FACTOR should be > 1.0) by \ 327squeezing in the frequency domain, \ 328then using the inverse DFT to get the time domain result.") 329 def stretch_sound_via_dft(factor, snd = false, chn = false) 330 factor = factor.to_f 331 n = framples(snd, chn) 332 n2 = (n / 2.0).floor 333 out_n = (n * factor).round 334 in_data = channel2vct(0, n, snd, chn) 335 out_data = make_vct(out_n) 336 fr = make_array(out_n, 0.0) 337 freq = (PI * 2) / n 338 n.times do |i| 339 d = edot_product(freq * Complex(0.0, 1.0) * i, in_data) 340 if i < n2 341 fr[i] = d 342 else 343 fr[i + (out_n - n - 1)] = d 344 end 345 end 346 freq = (PI * 2) / out_n 347 out_n.times do |i| 348 out_data[i] = (edot_product(freq * Complex(0.0, 1.0) * i, fr) / n).real 349 end 350 vct2channel(out_data, 0, out_n, snd, chn, false, 351 format("%s(%s", get_func_name, factor)) 352 end 353 354 # compute-uniform-circular-string 355 # 356 # this is a simplification of the underlying table-filling routine 357 # for "scanned synthesis". To watch the wave, open some sound (so 358 # Snd has some place to put the graph), turn off the time domain 359 # display (to give our graph all the window -- to do this in a much 360 # more elegant manner, see snd-motif.scm under scanned-synthesis). 361 def compute_uniform_circular_string(size, x0, x1, x2, mass, xspring, damp) 362 circle_vct_ref = lambda do |v, i| 363 if i < 0 364 v[i + size] 365 elsif i >= size 366 v[i - size] 367 else 368 v[i] 369 end 370 end 371 dm = damp / mass.to_f 372 km = xspring / mass.to_f 373 denom = 1.0 + dm 374 p1 = (2.0 + (dm - 2.0 * km)) / denom 375 p2 = km / denom 376 p3 = -1.0 / denom 377 size.times do |i| 378 x0[i] = p1 * x1[i] + 379 p2 * (circle_vct_ref.call(x1, i - 1) + 380 circle_vct_ref.call(x1, i + 1)) + 381 p3 * x2[i] 382 end 383 vct_fill!(x2, 0.0) 384 vct_add!(x2, x1) 385 vct_fill!(x1, 0.0) 386 vct_add!(x1, x0) 387 end 388 389 def testunif(mass, xspring, damp) 390 size = 128 391 x0 = make_vct(size) 392 x1 = make_vct(size) 393 x2 = make_vct(size) 394 12.times do |i| 395 x1[i + size / 4 - 6] = sin((TWO_PI * i) / 12.0) 396 end 397 1024.times do |i| 398 compute_uniform_circular_string(size, x0, x1, x2, mass, xspring, damp) 399 graph(x0, "string", 0, 1.0, -10.0, 10.0) 400 end 401 end 402 403 def test_scanned_synthesis(amp, dur, mass, xspring, damp) 404 size = 256 405 x0 = make_vct(size) 406 gx1 = make_vct(size) 407 gx2 = make_vct(size) 408 12.times do |i| 409 gx1[i + size / 4 - 6] = sin((TWO_PI * i) / 12.0) 410 end 411 gen1 = make_table_lookup(440.0, :wave, gx1) 412 gen2 = make_table_lookup(440.0, :wave, gx2) 413 x1 = gen1.data 414 x2 = gen2.data 415 recompute_samps = 30.0 416 k = 0.0 417 kincr = 1.0 / recompute_samps 418 data = make_vct!(dur) do |i| 419 if k >= 1.0 420 k = 0.0 421 compute_uniform_circular_string(size, x0, x1, x2, mass, xspring, damp) 422 else 423 k += kincr 424 end 425 g1 = table_lookup(gen1) 426 g2 = table_lookup(gen2) 427 g2 + k * (g1 - g2) 428 end 429 vct_scale!(data, amp / vct_peak(data)) 430 vct2channel(data, 0, dur) 431 end 432 433 # this is the more general form 434 def compute_string(size, x0, x1, x2, 435 masses, xsprings, esprings, damps, haptics) 436 circle_vct_ref = lambda do |v, i| 437 if i < 0 438 v[i + size] 439 elsif i >= size 440 v[i - size] 441 else 442 v[i] 443 end 444 end 445 size.times do |i| 446 dm = damps[i] / masses[i] 447 km = xsprings[i] / masses[i] 448 cm = esprings[i] / masses[i] 449 denom = 1.0 + dm + cm 450 p1 = (2.0 + (dm - 2.0 * km)) / denom 451 p2 = km / denom 452 p3 = -1.0 / denom 453 p4 = haptics / (masses[i] * denom) 454 x0[i] = p1 * x1[i] + 455 p2 * (circle_vct_ref.call(x1, i - 1) + 456 circle_vct_ref.call(x1, i + 1)) + 457 p3 * x2[i] + 458 p4 459 end 460 size.times do |i| 461 x2[i], x1[i] = x1[i], x0[i] 462 end 463 end 464 465 # "frequency division" -- an effect from sed_sed@my-dejanews.com 466 add_help(:freqdiv, 467 "freqdiv(n, snd=false, chn=false) \ 468Repeats each nth sample N times (clobbering the intermediate samples): \ 469freqdiv(8)") 470 def freqdiv(n, snd = false, chn = false) 471 div = 0 472 curval = 0.0 473 map_channel(lambda do |val| 474 curval = val if div.zero? 475 div += 1 476 div = 0 if div == n 477 curval 478 end, 0, false, snd, chn, false, 479 format("%s(%s", get_func_name, n)) 480 end 481 482 # "adaptive saturation" -- an effect from sed_sed@my-dejanews.com 483 # 484 # a more extreme effect is "saturation": 485 # (map-channel (lambda (val) 486 # (if (< (abs val) .1) val (if (>= val 0.0) 0.25 -0.25)))) 487 add_help(:adsat, 488 "adsat(size, beg=false, dur=false, snd=false, chn=false) \ 489Is an 'adaptive saturation' sound effect.") 490 def adsat(size, beg = false, dur = false, snd = false, chn = false) 491 mn = 0.0 492 mx = 0.0 493 n = 0 494 vals = make_vct(size) 495 map_channel(lambda do |val| 496 if n == size 497 size.times do |i| 498 if vals[i] >= 0.0 499 vals[i] = mx 500 else 501 vals[i] = mn 502 end 503 end 504 n = 0 505 mx = 0.0 506 mn = 0.0 507 vals 508 else 509 vals[n] = val 510 mx = val if val > mx 511 mn = val if val < mn 512 n += 1 513 false 514 end 515 end, beg, dur, snd, chn, false, 516 format("%s(%s, %s, %s", get_func_name, size, beg, dur)) 517 end 518 519 # spike 520 # 521 # makes sound more spikey -- sometimes a nice effect 522 add_help(:spike, 523 "spike(snd=false, chn=false) \ 524Multiplies successive samples together to make a sound more spikey.") 525 def spike(snd = false, chn = false) 526 x1 = x2 = 0.0 527 amp = maxamp(snd, chn) 528 map_channel(lambda do |x0| 529 res = (x0 / (amp * amp)) * x2.abs * x1.abs 530 x2, x1 = x1, x0 531 res 532 end, 0, false, snd, chn, false, "spike(") 533 end 534 535 # easily-fooled autocorrelation-based pitch tracker 536 add_help(:spot_freq, 537 "spot_freq(samp=0, snd=false, chn=false) \ 538Tries to determine the current pitch: spot_freq(left_sample)") 539 def spot_freq(samp = 0, snd = false, chn = false) 540 pow2 = (log(srate(snd) / 20.0) / log(2)).ceil 541 fftlen = (2 ** pow2).round 542 data = autocorrelate(channel2vct(samp, fftlen, snd, chn)) 543 cor_peak = vct_peak(data) 544 cor_peak2 = 2.0 * cor_peak 545 ret = 0.0 546 (1...fftlen - 2).each do |i| 547 if data[i] < data[i + 1] and data[i + 2] < data[i + 1] 548 logla = log10((cor_peak + data[i]) / cor_peak2) 549 logca = log10((cor_peak + data[i + 1]) / cor_peak2) 550 logra = log10((cor_peak + data[i + 2]) / cor_peak2) 551 offset = (0.5 * (logla - logra)) / (logla + logra + -2.0 * logca) 552 ret = srate(snd) / (2.0 * (i + 1 + offset)) 553 break 554 end 555 end 556 ret 557 end 558 # $graph_hook.add_hook!("examp-left-sample-hook") do |snd, chn, y0, y1| 559 # msg = format("(freq: %.3f)", spot_freq(left_sample(snd, chn))) 560 # status_report(msg, snd) 561 # end 562 # 563 # or 564 # 565 # $mouse_click_hook.add_hook!("examp-cursor-hook") do |snd, chn, 566 # button, state, 567 # x, y, axis| 568 # if axis == Time_graph 569 # status_report(format("(freq: %.3f)", spot_freq(cursor(snd, chn))), snd) 570 # end 571 # end 572 573 # chorus (doesn't always work and needs speedup) 574 class Flanger < Musgen 575 def initialize(time = 0.05, amount = 20.0, speed = 10.0) 576 super() 577 @time = time 578 @amount = amount 579 @speed = speed 580 @randind = make_rand_interp(:frequency, speed, :amplitude, amount) 581 @data = @randind.data 582 @length = @randind.length 583 len = random(3.0 * time * mus_srate()).floor 584 @flanger = make_delay(:size, len, :max_size, (len + amount + 1).to_i) 585 end 586 587 def inspect 588 format("%s.new(%s, %s, %s)", self.class, @time, @amount, @speed) 589 end 590 591 def to_s 592 format("#<%s time: %1.3f, amount: %1.3f, speed: %1.3f>", 593 self.class, @time, @amount, @speed) 594 end 595 596 def run_func(val1 = 0.0, val2 = 0.0) 597 flanger(val1) 598 end 599 600 def flanger(inval) 601 inval + delay(@flanger, inval, rand_interp(@randind)) 602 end 603 end 604 605 def make_flanger(time = 0.05, amount = 20.0, speed = 10.0) 606 Flanger.new(time, amount, speed) 607 end 608 609 def flanger?(obj) 610 obj.kind_of?(Flanger) 611 end 612 613 def flanger(gen, inval) 614 gen.flanger(inval) 615 end 616 617 add_help(:chorus, 618 "chorus(size=5) \ 619Tries to produce the chorus sound effect.") 620 def chorus(size = 5) 621 dlys = make_array(size) do 622 make_flanger 623 end 624 sum = 0.0 625 lambda do |inval| 626 dlys.each do |dly| 627 sum += dly.flanger(inval) 628 end 629 sum * 0.25 630 end 631 end 632 633 # chordalize (comb filters to make a chord using chordalize-amount 634 # and chordalize-base) 635 add_help(:chordalize, 636 "chordalize(amount=0.95, base=100, chord=[1.00, 0.75, 1.25]) \ 637Uses harmonically-related comb-filters to bring out a chord in a sound.") 638 def chordalize(amount = 0.95, base = 100, chord = [1.00, 0.75, 1.25]) 639 combs = chord.map do |interval| 640 make_comb(:scaler, amount, :size, (base * interval).round) 641 end 642 scaler = 0.5 / chord.length 643 lambda do |x| 644 val = 0.0 645 combs.each do |c| 646 val += comb(c, x) 647 end 648 scaler * val 649 end 650 end 651 652 # zero-phase, rotate-phase 653 # fft games (from the "phazor" package of Scott McNab) 654 add_help(:zero_phase, 655 "zero_phase(snd=false, chn=false) \ 656Calls fft, sets all phases to 0, and un-ffts.") 657 def zero_phase(snd = false, chn = false) 658 len = framples(snd, chn) 659 pow2 = (log(len) / log(2)).ceil 660 fftlen = (2 ** pow2).round 661 fftscale = 1.0 / fftlen 662 rl = channel2vct(0, fftlen, snd, chn) 663 old_pk = vct_peak(rl) 664 im = make_vct(fftlen) 665 fft(rl, im, 1) 666 rectangular2polar(rl, im) 667 vct_scale!(rl, fftscale) 668 vct_scale!(im, 0.0) 669 fft(rl, im, -1) 670 pk = vct_peak(rl) 671 vct2channel(rl.scale(old_pk / pk), 0, len, snd, chn, false, "zero_phase(") 672 end 673 674 # (set_)edit_list_proc_counter is defined in clm.rb 675 # It's necessary to produce a uniq method name. 676 add_help(:rotate_phase, 677 "rotate_phase(func, snd=false, chn=false) \ 678Calls fft, applies func to each phase, then un-ffts.") 679 def rotate_phase(func, snd = false, chn = false) 680 func_name = format("%s_%d", 681 get_func_name, 682 set_edit_list_proc_counter).intern 683 # Proc converted to Method (ie. normal function) for edit_list2function 684 func.to_method(func_name) 685 len = framples(snd, chn) 686 pow2 = (log(len) / log(2)).ceil 687 fftlen = (2 ** pow2).round 688 fftlen2 = (fftlen / 2).floor 689 fftscale = 1.0 / fftlen 690 rl = channel2vct(0, fftlen, snd, chn) 691 im = make_vct(fftlen) 692 old_pk = rl.peak 693 fft(rl, im, 1) 694 rectangular2polar(rl, im) 695 rl.scale!(fftscale) 696 im[0] = 0.0 697 j = fftlen - 1 698 (1...fftlen2).each do |i| 699 im[i] = snd_func(func_name, im[i]) 700 im[j] = -im[i] 701 j -= 1 702 end 703 polar2rectangular(rl, im) 704 fft(rl, im, -1) 705 pk = rl.peak 706 vct2channel(rl.scale(old_pk / pk), 0, len, snd, chn, false, 707 format("%s(Proc.new {|val_r| %s(val_r) }", 708 get_func_name, func_name)) 709 end 710 # rotate_phase(lambda {|x| 0.0 }) # is the same as (zero-phase) 711 # rotate_phase(lambda {|x| random(PI) }) # randomizes phases 712 # rotate_phase(lambda {|x| x }) # returns original 713 # rotate_phase(lambda {|x| -x }) # reverses original 714 # # (might want to write fftlen samps here) 715 716 # asymmetric FM (bes-i0 case) 717 class Asyfm < Musgen 718 def initialize(*args) 719 super() 720 frequency, ratio, r, index, phase = nil 721 optkey(args, binding, 722 [:frequency, $clm_default_frequency], 723 [:ratio, 1.0], 724 [:r, 1.0], 725 [:index, 1.0], 726 [:phase, 0.0]) 727 @frequency = frequency.to_f 728 @ratio = ratio.to_f 729 @r = r.to_f 730 @index = index.to_f 731 @freq = hz2radians(@frequency) 732 @phase = phase.to_f 733 end 734 attr_accessor :ratio, :r, :index 735 736 def inspect 737 format("%s.new(:frequency, %s, :ratio, %s, :r, %s, \ 738:index, %s, :freq, %s, :phase, %s)", 739 self.class, @frequency, @ratio, @r, @index, @freq, @phase) 740 end 741 742 def to_s 743 format("#<%s freq: %1.3f, ratio: %1.3f, r: %1.3f, \ 744index: %1.3f, freq: %1.3f, phase: %1.3f>", 745 self.class, @frequency, @ratio, @r, @index, @freq, @phase) 746 end 747 748 def run_func(val1 = 0.0, val2 = 0.0) 749 asyfm_J(val1) 750 end 751 752 def asyfm_J(input) 753 # It follows now asyfm-J in generators.scm, not dsp-asyfm-J in clm23.scm. 754 r1 = 1.0 / @r 755 one = ((@r > 1.0) or (@r < 0.0 and @r > -1.0)) ? -1.0 : 1.0 756 modphase = @ratio * @phase 757 result = exp(0.5 * @index * (@r - r1) * (one + cos(modphase))) * 758 cos(@phase + 0.5 * @index * (@r + r1) * sin(modphase)) 759 @phase = @phase + input + @freq 760 result 761 end 762 763 def asyfm_I(input) 764 r1 = 1.0 / @r 765 modphase = @ratio * @phase 766 result = exp(0.5 * @index * (@r + r1) * (cos(modphase) - 1.0)) - 767 cos(@phase + 0.5 * @index * (@r - r1) * sin(modphase)) 768 @phase = @phase + input + @freq 769 result 770 end 771 end 772 773 def make_asyfm(*args) 774 Asyfm.new(*args) 775 end 776 777 def asyfm?(obj) 778 obj.kind_of?(Asyfm) 779 end 780 781 def asyfm_J(gen, input) 782 gen.asyfm_J(input) 783 end 784 785 def asyfm_I(gen, input) 786 gen.asyfm_I(input) 787 end 788 789 # cosine-summation (a simpler version of sine-summation) 790 # 791 # from Andrews, Askey, Roy "Special Functions" 5.1.16 792 add_help(:cosine_summation, 793 "cosine_summation(gen, r) \ 794Is a variant of the CLM sine-summation generator; \ 795R controls successive sinusoid amplitudes.") 796 def cosine_summation(gen, r) 797 rr = r * r 798 rrp1 = 1.0 + rr 799 rrm1 = 1.0 - rr 800 r2 = 2.0 * r 801 ((rrm1 / (rrp1 - r2 * oscil(gen))) - 1.0) * ((1.0 - r) / r2) 802 end 803 alias make_cosine_summation make_oscil 804 805 # kosine-summation 806 # 807 # from Askey "Ramanujan and Hypergeometric Series" in Berndt and 808 # Rankin "Ramanujan: Essays and Surveys" 809 # 810 # this gives a sum of cosines of decreasing amp where the "k" 811 # parameter determines the "index" (in FM nomenclature) -- higher k 812 # = more cosines; the actual amount of the nth cos involves 813 # hypergeometric series (looks like r^n/n! (~=e^n?) with a million 814 # other terms). 815 add_help(:kosine_summation, 816 "kosine_summation(gen, r, k) \ 817Is a variant of sum-of-cosines; \ 818R controls successive sinusoid amplitude; \ 819K controls how many sinusoids are produced.") 820 def kosine_summation(gen, r, k) 821 r2 = r * r 822 ((1.0 + r2) - 2 * r * oscil(gen)) ** -k * ((1.0 + r2) - 2 * r) ** k 823 end 824 alias make_kosine_summation make_oscil 825 826 # legendre, fejer 827 def fejer_sum(angle, n) 828 if angle.zero? 829 1.0 830 else 831 val = sin(0.5 * (n + 1) * angle) / (2.0 * sin(0.5 * angle)) 832 2.0 * ((val * val) / (n + 1)) 833 end 834 end 835 836 def legendre_sum(angle, n) 837 val = sin(angle * (n + 0.5)) / sin(0.5 * angle) 838 val * val 839 end 840 841 # variations on sum-of-cosines 842 # from "Trigonometric Delights" by Eli Maor 843 def sum_of_n_sines(angle, n) 844 a2 = angle * 0.5 845 den = sin(a2) 846 if den.zero? 847 0.0 848 else 849 (sin(n * a2) * sin((n + 1) * a2)) / den 850 end 851 end 852 853 def sum_of_n_odd_sines(angle, n) 854 angle = angle.to_f 855 n = n.to_f 856 den = sin(angle) 857 na = sin(n * angle) 858 if den.zero? 859 0.0 860 else 861 (na * na) / den 862 end 863 end 864 865 def sum_of_n_odd_cosines(angle, n) 866 angle = angle.to_f 867 n = n.to_f 868 den = 2.0 * sin(angle) 869 if den.zero? 870 n 871 else 872 sin(2.0 * n * angle) / den 873 end 874 end 875 876 # x = current phase, a = amp (more or less), N = 1..10 or 877 # thereabouts, fi = phase increment Alexander Kritov suggests 878 # time-varying "a" is good (this is a translation of his code) from 879 # Stilson/Smith apparently -- was named "Discrete Summation Formula" 880 # which doesn't convey anything to me 881 def band_limited_sawtooth(x, a, n, fi) 882 s4 = 1.0 + -2.0 * a * cos(x) + a * a 883 if s4.zero? 884 0.0 885 else 886 s1 = a ** (n - 1.0) * sin((n - 1.0) * x + fi) 887 s2 = a ** n * sin(n * x + fi) 888 s3 = a * sin(x + fi) 889 (sin(fi) + -s3 + -s2 + s1) / s4 890 end 891 end 892 893 # brighten-slightly 894 add_help(:brighten_slightly, 895 "brighten_slightly(amount, snd=false, chn=false) \ 896Is a form of contrast-enhancement (AMOUNT between ca 0.1 and 1.0).") 897 def brighten_slightly(amount, snd = false, chn = false) 898 mx = maxamp 899 brt = (TWO_PI * amount) / mx 900 map_channel(lambda do |y| 901 mx * sin(y * brt) 902 end, 0, false, snd, chn, false, 903 format("%s(%s", get_func_name, amount)) 904 end 905 906 add_help(:brighten_slightly_1, 907 "brighten_slightly_1(coeffs, snd=false, chn=false) \ 908Is a form of contrast-enhancement: brighten_slightly-1([1, 0.5, 3, 1])") 909 def brighten_slightly_1(coeffs, snd = false, chn = false) 910 pcoeffs = partials2polynomial(coeffs) 911 mx = maxamp(snd, chn) 912 map_channel(lambda do |y| 913 mx * polynomial(pcoeffs, y / mx) 914 end, 0, false, snd, chn, false, 915 format("%s(%s", get_func_name, coeffs)) 916 end 917 918 # FIR filters 919 920 # Snd's (very simple) spectrum->coefficients procedure is: 921 add_help(:spectrum2coeffs, 922 "spectrum2coeffs(order, spectr) \ 923Returns FIR filter coefficients given the filter ORDER \ 924and desired spectral envelope (a vct).") 925 def spectrum2coeffs(order, spectr) 926 coeffs = make_vct(order) 927 n = order 928 m = ((n + 1) / 2.0).floor 929 am = 0.5 * (n + 1) 930 q = (PI * 2.0) / n 931 jj = n - 1 932 m.times do |j| 933 xt = 0.5 * spectr[0] 934 (1...m).each do |i| 935 xt = xt + spectr[i] * cos(q * i * (am - j - 1)) 936 end 937 coeff = 2.0 * (xt / n) 938 coeffs[j] = coeff 939 coeffs[jj] = coeff 940 jj -= 1 941 end 942 coeffs 943 end 944 945 add_help(:fltit_1, 946 "fltit_1(order, spectrum) \ 947Creates an FIR filter from SPECTRUM and ORDER and \ 948returns a closure that calls it: \ 949map_channel(fltit_1(10, vct(0, 1.0, 0, 0, 0, 0, 0, 0, 1.0, 0)))") 950 def fltit_1(order, spectr) 951 flt = make_fir_filter(order, spectrum2coeffs(order, spectr)) 952 lambda do |x| 953 fir_filter(flt, x) 954 end 955 end 956 957 # Hilbert transform 958 add_help(:make_hilbert_transform, 959 "make_hilbert_transform(len=30) \ 960Makes a Hilbert transform filter.") 961 def make_hilbert_transform(len = 30) 962 arrlen = len * 2 + 1 963 arr = make_vct(arrlen) 964 (-len...len).each do |i| 965 k = i + len 966 denom = PI * i 967 num = 1.0 - cos(denom) 968 if num.zero? or i.zero? 969 arr[k] = 0.0 970 else 971 arr[k] = (num / denom) * (0.54 + 0.46 * cos(denom / len)) 972 end 973 end 974 make_fir_filter(arrlen, arr) 975 end 976 977 add_help(:hilbert_transform, 978 "hilbert_transform(f, in) \ 979Is the generator corresponding to make_hilbert_transform.") 980 alias hilbert_transform fir_filter 981 982 # highpass filter 983 add_help(:make_highpass, 984 "make_highpass(fc, len=30) \ 985Makes an FIR highpass filter.") 986 def make_highpass(fc, len = 30) 987 fc = fc.to_f 988 arrlen = len * 2 + 1 989 arr = make_vct(arrlen) 990 (-len...len).each do |i| 991 k = i + len 992 denom = PI * i 993 num = -sin(fc * i) 994 if i.zero? 995 arr[k] = 1.0 - fc / PI 996 else 997 arr[k] = (num / denom) * (0.54 + 0.46 * cos((PI * i) / len)) 998 end 999 end 1000 make_fir_filter(arrlen, arr) 1001 end 1002 1003 add_help(:highpass, 1004 "highpass(f, in) \ 1005Is the generator corresponding to make_highpass.") 1006 alias highpass fir_filter 1007 1008 # lowpass filter 1009 add_help(:make_lowpass, 1010 "make_lowpass(fc, len=30) \ 1011Makes an FIR lowpass filter.") 1012 def make_lowpass(fc, len = 30) 1013 fc = fc.to_f 1014 arrlen = len * 2 + 1 1015 arr = make_vct(arrlen) 1016 (-len...len).each do |i| 1017 k = i + len 1018 denom = PI * i 1019 num = sin(fc * i) 1020 if i.zero? 1021 arr[k] = fc / PI 1022 else 1023 arr[k] = (num / denom) * (0.54 + 0.46 * cos((PI * i) / len)) 1024 end 1025 end 1026 make_fir_filter(arrlen, arr) 1027 end 1028 1029 add_help(:lowpass, 1030 "lowpass(f, in) \ 1031Is the generator corresponding to make_lowpass.") 1032 alias lowpass fir_filter 1033 1034 # bandpass filter 1035 add_help(:make_bandpass, 1036 "make_bandpass(flo, fhi, len=30) \ 1037Makes an FIR bandpass filter.") 1038 def make_bandpass(flo, fhi, len = 30) 1039 flo = flo.to_f 1040 fhi = fhi.to_f 1041 arrlen = len * 2 + 1 1042 arr = make_vct(arrlen) 1043 (-len...len).each do |i| 1044 k = i + len 1045 denom = PI * i 1046 num = sin(fhi * i) - sin(flo * i) 1047 if i.zero? 1048 arr[k] = (fhi - flo) / PI 1049 else 1050 arr[k] = (num / denom) * (0.54 + 0.46 * cos((PI * i) / len)) 1051 end 1052 end 1053 make_fir_filter(arrlen, arr) 1054 end 1055 1056 add_help(:bandpass, 1057 "bandpass(f, in) \ 1058Is the generator corresponding to make_bandpass.") 1059 alias bandpass fir_filter 1060 1061 # bandstop filter 1062 add_help(:make_bandstop, 1063 "make_bandstop(flo, fhi, len=30) \ 1064Makes an FIR bandstop (notch) filter.") 1065 def make_bandstop(flo, fhi, len = 30) 1066 flo = flo.to_f 1067 fhi = fhi.to_f 1068 arrlen = len * 2 + 1 1069 arr = make_vct(arrlen) 1070 (-len...len).each do |i| 1071 k = i + len 1072 denom = PI * i 1073 num = sin(flo * i) - sin(fhi * i) 1074 if i.zero? 1075 arr[k] = 1.0 - (fhi - flo) / PI 1076 else 1077 arr[k] = (num / denom) * (0.54 + 0.46 * cos((PI * i) / len)) 1078 end 1079 end 1080 make_fir_filter(arrlen, arr) 1081 end 1082 1083 add_help(:bandstop, 1084 "bandstop(f, in) \ 1085Is the generator corresponding to make_bandstop.") 1086 alias bandstop fir_filter 1087 1088 # differentiator 1089 add_help(:make_differentiator, 1090 "make_differentiator(len=30) \ 1091Makes an FIR differentiator (highpass) filter.") 1092 def make_differentiator(len = 30) 1093 arrlen = len * 2 + 1 1094 arr = make_vct(arrlen) 1095 (-len...len).each do |i| 1096 k = i + len 1097 if i.zero? 1098 arr[k] = 0.0 1099 else 1100 arr[k] = (cos(PI * i) / i - sin(PI * i) / (PI * i * i)) * 1101 (0.54 + 0.46 * cos((PI * i) / len)) 1102 end 1103 end 1104 make_fir_filter(arrlen, arr) 1105 end 1106 1107 add_help(:differentiator, 1108 "differentiator(f, in) \ 1109Is the generator corresponding to make_differentiator.") 1110 alias differentiator fir_filter 1111 1112 # IIR filters 1113 # 1114 # Butterworth filters (see also further below -- make-butter-lp et al) 1115 # 1116 # translated from CLM butterworth.cl: 1117 # 1118 # Sam Heisz, January 1998 1119 # inspired by some unit generators written for Csound by Paris 1120 # Smaragdis who based his work on formulas from Charles Dodge, 1121 # Computer music: synthesis, composition, and performance. 1122 1123 add_help(:butter, 1124 "butter(b, sig) \ 1125Is the generator side for the various make_butter procedure.") 1126 alias butter filter 1127 add_help(:make_butter_high_pass, 1128 "make_butter_high_pass(freq) \ 1129Makes a Butterworth filter with high pass cutoff at FREQ.") 1130 def make_butter_high_pass(freq) 1131 r = tan(PI * freq / srate()) 1132 r2 = r * r 1133 c1 = 1.0 / (1.0 + r * sqrt(2.0) + r2) 1134 c2 = -2.0 * c1 1135 c3 = c1 1136 c4 = 2.0 * (r2 - 1.0) * c1 1137 c5 = ((1.0 - r * sqrt(2.0)) + r2) * c1 1138 make_filter(3, vct(c1, c2, c3), vct(0.0, c4, c5)) 1139 end 1140 1141 add_help(:make_butter_low_pass, 1142 "make_butter_low_pass(freq) \ 1143Makes a Butterworth filter with high pass cutoff at FREQ. \ 1144The result can be used directly: \ 1145filter_sound(make_butter_low_pass(500.0)), or via the 'butter' generator.") 1146 def make_butter_low_pass(freq) 1147 r = 1.0 / tan(PI * freq / srate()) 1148 r2 = r * r 1149 c1 = 1.0 / (1.0 + r * sqrt(2.0) + r2) 1150 c2 = 2.0 * c1 1151 c3 = c1 1152 c4 = 2.0 * (1.0 - r2) * c1 1153 c5 = ((1.0 - r * sqrt(2.0)) + r2) * c1 1154 make_filter(3, vct(c1, c2, c3), vct(0.0, c4, c5)) 1155 end 1156 1157 add_help(:make_butter_band_pass, 1158 "make_butter_band_pass(freq, band) \ 1159Makes a bandpass Butterworth filter with low edge at FREQ and width BAND.") 1160 def make_butter_band_pass(freq, band) 1161 d = 2.0 * cos(2.0 * PI * freq / srate()) 1162 c = 1.0 / tan(PI * band / srate()) 1163 c1 = 1.0 / (1.0 + c) 1164 c2 = 0.0 1165 c3 = -c1 1166 c4 = -c * d * c1 1167 c5 = (c - 1.0) * c1 1168 make_filter(3, vct(c1, c2, c3), vct(0.0, c4, c5)) 1169 end 1170 1171 add_help(:make_butter_band_reject, 1172 "make_butter_band_reject(freq, band) \ 1173Makes a band-reject Butterworth filter with low edge at FREQ and width BAND.") 1174 def make_butter_band_reject(freq, band) 1175 d = 2.0 * cos(2.0 * PI * freq / srate()) 1176 c = tan(PI * band / srate()) 1177 c1 = 1.0 / (1.0 + c) 1178 c2 = -d * c1 1179 c3 = c1 1180 c4 = c2 1181 c5 = (1.0 - c) * c1 1182 make_filter(3, vct(c1, c2, c3), vct(0.0, c4, c5)) 1183 end 1184 1185 # from "DSP Filter Cookbook" by Lane et al, Prompt Pubs, 2001 1186 # 1187 # use with the filter generator 1188 # (define gen (make-iir-high-pass-2 1000)) 1189 # (filter gen 1.0) 1190 # etc 1191 add_help(:make_biquad, 1192 "make_biquad(a0, a1, a2, b1, b2) \ 1193Returns a biquad filter (use with the CLM filter gen).") 1194 def make_biquad(a0, a1, a2, b1, b2) 1195 make_filter(3, vct(a0, a1, a2), vct(0.0, b1, b2)) 1196 end 1197 1198 # din=(sqrt 2.0) for example (suggested range 0.2...10) 1199 def make_iir_low_pass_2(fc, din = false) 1200 fc = fc.to_f 1201 theta = (TWO_PI * fc) / mus_srate() 1202 d = (din or sqrt(2.0)) 1203 beta = 0.5 * ((1.0 - (d / 2.0) * sin(theta)) / 1204 (1.0 + (d / 2.0) * sin(theta))) 1205 gamma = (0.5 + beta) * cos(theta) 1206 alpha = 0.5 * (0.5 + beta + -gamma) 1207 make_filter(3, 1208 vct(alpha, 2.0 * alpha, alpha), 1209 vct(0.0, -2.0 * gamma, 2.0 * beta)) 1210 end 1211 1212 def make_iir_high_pass_2(fc, din = false) 1213 fc = fc.to_f 1214 theta = (TWO_PI * fc) / mus_srate() 1215 d = (din or sqrt(2.0)) 1216 beta = 0.5 * ((1.0 - (d / 2.0) * sin(theta)) / 1217 (1.0 + (d / 2.0) * sin(theta))) 1218 gamma = (0.5 + beta) * cos(theta) 1219 alpha = 0.5 * (0.5 + beta + gamma) 1220 make_filter(3, 1221 vct(alpha, -2.0 * alpha, alpha), 1222 vct(0.0, -2.0 * gamma, 2.0 * beta)) 1223 end 1224 1225 def make_iir_band_pass_2(f1, f2) 1226 f1 = f1.to_f 1227 f2 = f2.to_f 1228 theta = (TWO_PI * sqrt(f1 * f2)) / mus_srate() 1229 q = sqrt(f1 * f2) / (f2 - f1) 1230 t2 = tan(theta / (2 * q)) 1231 beta = 0.5 * ((1.0 - t2) / (1.0 + t2)) 1232 gamma = (0.5 + beta) * cos(theta) 1233 alpha = 0.5 - beta 1234 make_filter(3, 1235 vct(alpha, 0.0, -alpha), 1236 vct(0.0, -2.0 * gamma, 2.0 * beta)) 1237 end 1238 1239 def make_iir_band_stop_2(f1, f2) 1240 f1 = f1.to_f 1241 f2 = f2.to_f 1242 theta = (TWO_PI * sqrt(f1 * f2)) / mus_srate() 1243 q = sqrt(f1 * f2) / (f2 - f1) 1244 t2 = tan(theta / (2 * q)) 1245 beta = 0.5 * ((1.0 - t2) / (1.0 + t2)) 1246 gamma = (0.5 + beta) * cos(theta) 1247 alpha = 0.5 + beta 1248 make_filter(3, 1249 vct(alpha, -2.0 * gamma, alpha), 1250 vct(0.0, -2.0 * gamma, 2.0 * beta)) 1251 end 1252 1253 def make_eliminate_hum(hum_freq = 60.0, hum_harmonics = 5, bandwidth = 10) 1254 b2 = 0.5 * bandwidth 1255 make_array(hum_harmonics) do |i| 1256 center = (i + 1.0) * hum_freq 1257 make_iir_band_stop_2(center - b2, center + b2) 1258 end 1259 end 1260 1261 def eliminate_hum(gens, x0) 1262 val = x0 1263 gens.each do |gen| 1264 val = filter(gen, val) 1265 end 1266 val 1267 end 1268 1269 # bandpass, m is gain at center of peak 1270 # use map-channel with this one (not clm-channel or filter) 1271 def make_peaking_2(f1, f2, m) 1272 f1 = f1.to_f 1273 f2 = f2.to_f 1274 theta = (TWO_PI * sqrt(f1 * f2)) / mus_srate() 1275 q = sqrt(f1 * f2) / (f2 - f1) 1276 t2 = (4.0 / (m + 1.0)) * tan(theta / (2 * q)) 1277 beta = 0.5 * ((1.0 - t2) / (1.0 + t2)) 1278 gamma = (0.5 + beta) * cos(theta) 1279 alpha = 0.5 - beta 1280 flt = make_filter(3, 1281 vct(alpha, 0.0, -alpha), 1282 vct(0.0, -2.0 * gamma, 2.0 * beta)) 1283 lambda do |x| 1284 x + (m - 1.0) * filter(flt, x) 1285 end 1286 end 1287 1288 # convert cascade coeffs to canonical form 1289 # from Orfanidis "Introduction to Signal Processing 1290 def c2c_conv(m, h, l, x, y) 1291 (l + m).times do |i| 1292 y[i] = 0.0 1293 ([0, i - (1 + l)].max..[i, m].min).each do |j| 1294 y[i] = y[i] + h[j] * x[i - j] 1295 end 1296 end 1297 end 1298 1299 add_help(:cascade2canonical, 1300 "cascade2canonical(a) \ 1301Converts a list of cascade coeffs (vcts with 3 entries) to canonical form.") 1302 def cascade2canonical(a) 1303 k = a.length 1304 d = make_vct(2 * k + 1) 1305 a1 = make_vct(2 * k + 1) 1306 a1[0] = 1.0 1307 k.times do |i| 1308 c2c_conv(2, a[i], 2 * i + 1, a1, d) 1309 (2 * i + 3).times do |j| 1310 a1[j] = d[j] 1311 end 1312 end 1313 a1 1314 end 1315 1316 # order is M*2, fc is cutoff freq (Hz) 1317 add_help(:make_butter_lp, 1318 "make_butter_lp(m, fc) \ 1319Returns a butterworth low-pass filter; \ 1320its order is M * 2, FC is the cutoff frequency in Hz.") 1321 def make_butter_lp(m, fc) 1322 fc = fc.to_f 1323 xcoeffs = make_array(m) 1324 ycoeffs = make_array(m) 1325 theta = (TWO_PI * fc) / mus_srate() 1326 st = sin(theta) 1327 ct = cos(theta) 1328 m.times do |k| 1329 d = 2.0 * sin((PI * (2.0 * (k + 1.0) - 1.0)) / (4.0 * m)) 1330 beta = 0.5 * ((1.0 - 0.5 * d * st) / (1.0 + 0.5 * d * st)) 1331 gamma = ct * (0.5 + beta) 1332 alpha = 0.25 * (0.5 + beta + -gamma) 1333 xcoeffs[k] = vct(2.0 * alpha, 4.0 * alpha, 2.0 * alpha) 1334 ycoeffs[k] = vct(1.0, -2.0 * gamma, 2.0 * beta) 1335 end 1336 make_filter(2 * m + 1, 1337 cascade2canonical(xcoeffs), cascade2canonical(ycoeffs)) 1338 end 1339 1340 # order is M*2, fc is cutoff freq (Hz) 1341 add_help(:make_butter_hp, 1342 "make_butter_hp(m, fc) \ 1343Returns a butterworth high-pass filter; \ 1344its order is M * 2, FC is the cutoff frequency in Hz.") 1345 def make_butter_hp(m, fc) 1346 fc = fc.to_f 1347 xcoeffs = make_array(m) 1348 ycoeffs = make_array(m) 1349 theta = (TWO_PI * fc) / mus_srate() 1350 st = sin(theta) 1351 ct = cos(theta) 1352 m.times do |k| 1353 d = 2.0 * sin((PI * (2.0 * (k + 1.0) - 1.0)) / (4.0 * m)) 1354 beta = 0.5 * ((1.0 - 0.5 * d * st) / (1.0 + 0.5 * d * st)) 1355 gamma = ct * (0.5 + beta) 1356 alpha = 0.25 * (0.5 + beta + gamma) 1357 xcoeffs[k] = vct(2.0 * alpha, -4 * alpha, 2.0 * alpha) 1358 ycoeffs[k] = vct(1.0, -2.0 * gamma, 2.0 * beta) 1359 end 1360 make_filter(2 * m + 1, 1361 cascade2canonical(xcoeffs), cascade2canonical(ycoeffs)) 1362 end 1363 1364 # order is M*2, f1 and f2 are band edge freqs (Hz) 1365 add_help(:make_butter_bp, 1366 "make_butter_bp(m, f1, f2) \ 1367Returns a butterworth band-pass filter; \ 1368its order is M * 2, F1 and F2 are the band edge frequencies in Hz.") 1369 def make_butter_bp(m, f1, f2) 1370 f1 = f1.to_f 1371 f2 = f2.to_f 1372 xcoeffs = make_array(m) 1373 ycoeffs = make_array(m) 1374 f0 = sqrt(f1 * f2) 1375 q = f0 / (f2 - f1) 1376 theta0 = (TWO_PI * f0) / mus_srate() 1377 de = (2.0 * tan(theta0 / (2.0 * q))) / sin(theta0) 1378 de2 = de / 2.0 1379 tn0 = tan(theta0 * 0.5) 1380 k = j = 1 1381 m.times do |i| 1382 dk = 2.0 * sin((PI * (2.0 * k - 1.0)) / (2.0 * m)) 1383 ak = (1.0 + de2 * de2) / (dk * de2) 1384 dk1 = sqrt((de * dk) / (ak + sqrt(ak * ak - 1.0))) 1385 bk = de2 * (dk / dk1) 1386 wk = (bk + sqrt(bk * bk - 1.0)).real 1387 thetajk = ((j == 1) ? (2.0 * atan(tn0 / wk)) : (2.0 * atan(tn0 * wk))) 1388 betajk = 0.5 * ((1.0 - 0.5 * dk1 * sin(thetajk)) / 1389 (1.0 + 0.5 * dk1 * sin(thetajk))) 1390 gammajk = (0.5 + betajk) * cos(thetajk) 1391 wk2 = (wk - 1.0 / wk) / dk1 1392 alphajk = 0.5 * (0.5 - betajk) * sqrt(1.0 + wk2 * wk2) 1393 xcoeffs[i] = vct(2.0 * alphajk, 0.0, -2.0 * alphajk) 1394 ycoeffs[i] = vct(1.0, -2.0 * gammajk, 2.0 * betajk) 1395 if j == 1 1396 j = 2 1397 else 1398 k += 1 1399 j = 1 1400 end 1401 end 1402 make_filter(2 * m + 1, 1403 cascade2canonical(xcoeffs), cascade2canonical(ycoeffs)) 1404 end 1405 1406 # order is M*2, f1 and f2 are band edge freqs (Hz) 1407 add_help(:make_butter_bs, 1408 "make_butter_bs(m, f1, f2) \ 1409Returns a butterworth band-stop filter; \ 1410its order is M * 2, F1 and F2 are the band edge frequencies in Hz.") 1411 def make_butter_bs(m, f1, f2) 1412 f1 = f1.to_f 1413 f2 = f2.to_f 1414 xcoeffs = make_array(m) 1415 ycoeffs = make_array(m) 1416 f0 = sqrt(f1 * f2) 1417 q = f0 / (f2 - f1) 1418 theta0 = (TWO_PI * f0) / mus_srate() 1419 de = (2.0 * tan(theta0 / (2.0 * q))) / sin(theta0) 1420 de2 = de / 2.0 1421 ct = cos(theta0) 1422 tn0 = tan(theta0 * 0.5) 1423 k = j = 1 1424 m.times do |i| 1425 dk = 2.0 * sin((PI * (2.0 * k - 1.0)) / (2.0 * m)) 1426 ak = (1.0 + de2 * de2) / (dk * de2) 1427 dk1 = sqrt((de * dk) / (ak + sqrt(ak * ak - 1.0))) 1428 bk = de2 * (dk / dk1) 1429 wk = (bk + sqrt(bk * bk - 1.0)).real 1430 thetajk = ((j == 1) ? (2.0 * atan(tn0 / wk)) : (2.0 * atan(tn0 * wk))) 1431 betajk = 0.5 * ((1.0 - 0.5 * dk1 * sin(thetajk)) / 1432 (1.0 + 0.5 * dk1 * sin(thetajk))) 1433 gammajk = (0.5 + betajk) * cos(thetajk) 1434 alphajk = 0.5 * (0.5 + betajk) * ((1.0 - cos(thetajk)) / (1.0 - ct)) 1435 xcoeffs[i] = vct(2.0 * alphajk, -4.0 * ct * alphajk, 2.0 * alphajk) 1436 ycoeffs[i] = vct(1.0, -2.0 * gammajk, 2.0 * betajk) 1437 if j == 1 1438 j = 2 1439 else 1440 k += 1 1441 j = 1 1442 end 1443 end 1444 make_filter(2 * m + 1, 1445 cascade2canonical(xcoeffs), cascade2canonical(ycoeffs)) 1446 end 1447 1448 # notch filters 1449 def make_notch_frequency_response(cur_srate, freqs, notch_width = 2) 1450 cur_srate = cur_srate.to_f 1451 notch_width = notch_width.to_f 1452 freq_response = [0.0, 1.0] 1453 freqs.each do |f| 1454 # left upper y hz 1455 freq_response.push((2.0 * (f - notch_width)) / cur_srate) 1456 # left upper y resp 1457 freq_response.push(1.0) 1458 # left bottom y hz 1459 freq_response.push((2.0 * (f - notch_width / 2.0)) / cur_srate) 1460 # left bottom y resp 1461 freq_response.push(0.0) 1462 # right bottom y hz 1463 freq_response.push((2.0 * (f + notch_width / 2.0)) / cur_srate) 1464 # right bottom y resp 1465 freq_response.push(0.0) 1466 # right upper y hz 1467 freq_response.push((2.0 * (f + notch_width)) / cur_srate) 1468 # right upper y resp 1469 freq_response.push(1.0) 1470 end 1471 freq_response.push(1.0, 1.0) 1472 end 1473 1474 add_help(:notch_channel, 1475 "notch_channel(freqs, order=false, beg=false, dur=false, \ 1476snd=false, chn=false, edpos=false, trunc=true, notch_width=2) \ 1477Returns notch filter removing freqs.") 1478 def notch_channel(freqs, 1479 filter_order = false, 1480 beg = false, 1481 dur = false, 1482 snd = false, 1483 chn = false, 1484 edpos = false, 1485 truncate = true, 1486 notch_width = 2) 1487 sr = srate(snd).to_f 1488 lm = [framples(snd, chn), 2 ** (log(sr / notch_width) / log(2.0)).floor].min 1489 filter_channel(make_notch_frequency_response(sr, freqs, notch_width), 1490 (filter_order or lm), 1491 beg, dur, snd, chn, edpos, truncate, 1492 format("%s(%p, %s, %s, %s", 1493 get_func_name, freqs, filter_order, beg, dur)) 1494 end 1495 1496 add_help(:notch_sound, 1497 "notch_sound(freqs, order=false, \ 1498snd=false, chn=false, notch_width=2) \ 1499Returns notch filter removing freqs.") 1500 def notch_sound(freqs, filter_order = false, 1501 snd = false, chn = false, notch_width = 2) 1502 sr = srate(snd).to_f 1503 lm = [framples(snd, chn), 2 ** (log(sr / notch_width) / log(2.0)).floor].min 1504 filter_sound(make_notch_frequency_response(sr, freqs, notch_width), 1505 (filter_order or lm), 1506 snd, chn, false, 1507 format("%s(%p, %s, 0, false", 1508 get_func_name, freqs, filter_order)) 1509 end 1510 1511 add_help(:notch_selection, 1512 "notch_selection(freqs, order=false, notch_width=2) \ 1513Returns notch filter removing freqs.") 1514 def notch_selection(freqs, filter_order = false, 1515 snd = false, chn = false, notch_width = 2) 1516 if selection? 1517 sr = selection_srate.to_f 1518 fr = selection_framples() 1519 lm = [fr, 2 ** (log(sr / notch_width) / log(2.0)).floor].min 1520 filter_selection(make_notch_frequency_response(sr, freqs, notch_width), 1521 (filter_order or lm)) 1522 end 1523 end 1524 1525 # fractional Fourier Transform, z transform 1526 # 1527 # translated from the fxt package of Joerg Arndt 1528 add_help(:fractional_fourier_transform, 1529 "fractional_fourier_transform(real, imaginary, n, angle) \ 1530Performs a fractional Fourier transform on data; \ 1531if angle=1.0, you get a normal Fourier transform.") 1532 def fractional_fourier_transform(fr, fi, n, v) 1533 hr = make_vct(n) 1534 hi = make_vct(n) 1535 ph0 = (v * TWO_PI) / n 1536 n.times do |w| 1537 sr = 0.0 1538 si = 0.0 1539 n.times do |k| 1540 phase = ph0 * k * w 1541 c = cos(phase) 1542 s = sin(phase) 1543 x = fr[k] 1544 y = fi[k] 1545 r = x * c - y * s 1546 i = y * c + x * s 1547 sr += r 1548 si += i 1549 hr[w] = sr 1550 hi[w] = si 1551 end 1552 end 1553 [hr, hi] 1554 end 1555 1556 # using vector to allow complex sums (z=e^2*pi*i/n -> fourier transform) 1557 # z_transform(data, n, exp(Complex(0.0, (2.0 / n) * PI))) 1558 add_help(:z_transform, 1559 "z_transform(data, n, z) \ 1560Performs a Z transform on data; \ 1561if z=e^2*pi*j/n you get a Fourier transform; \ 1562complex results in returned vector.") 1563 def z_transform(f, n, z) 1564 make_array(n) do |w| 1565 sum = 0.0 1566 t = 1.0 1567 m = z ** w 1568 n.times do |k| 1569 sum = sum + f[k] * t 1570 t *= m 1571 end 1572 sum 1573 end 1574 end 1575 1576 # slow Hartley transform 1577 # 1578 # taken from Perry Cook's SignalProcessor.m (the slow version of the 1579 # Hartley transform) 1580 add_help(:dht, 1581 "dht(data) \ 1582Returns the Hartley transform of DATA.") 1583 def dht(data) 1584 len = data.length 1585 arr = make_vct(len) 1586 w = TWO_PI / len 1587 len.times do |i| 1588 data.each_with_index do |val, j| 1589 arr[i] = arr[i] + val * (cos(i * j * w) + sin(i * j * w)) 1590 end 1591 end 1592 arr 1593 end 1594 1595 add_help(:find_sine, 1596 "find_sine(freq, beg, dur, snd=false) \ 1597Returns the amplitude and initial-phase (for sin) at FREQ between BEG and DUR.") 1598 def find_sine(freq, beg, dur, snd = false) 1599 incr = (TWO_PI * freq) / srate(snd) 1600 sw = 0.0 1601 cw = 0.0 1602 reader = make_sampler(beg, snd) 1603 dur.times do |i| 1604 samp = next_sample(reader) 1605 inc = i * incr 1606 sw = sw + samp * sin(inc) 1607 cw = cw + samp * cos(inc) 1608 end 1609 [2.0 * (sqrt(sw * sw + cw * cw) / dur), atan2(cw, sw)] 1610 end 1611 1612# this is a faster version of find-sine using the "Goertzel algorithm" 1613# taken from R Lyons "Understanding DSP" p 529 1614# it returns the same result as find_sine above if you take (* 2 (/ 1615# (goertzel...) dur)) -- see snd-test.rb examples 1616 1617 add_help(:goertzel, 1618 "goertzel(freq, beg=0, dur=false, snd=false) \ 1619Returns the amplitude of the FREQ spectral component.") 1620 def goertzel(freq, beg = 0, dur = false, snd = false) 1621 y0 = 0.0 1622 y1 = 0.0 1623 y2 = 0.0 1624 rfreq = (TWO_PI * freq) / srate(snd) 1625 cs = 2.0 * cos(rfreq) 1626 scan_channel(lambda do |y| 1627 y2, y1 = y1, y0 1628 y0 = (y1 * cs - y2) + y 1629 false 1630 end, beg, (dur or framples(snd)), snd) 1631 (y0 - y1 * exp(Complex(0.0, -rfreq))).abs 1632 end 1633 1634 add_help(:make_spencer_filter, 1635 "make_spencer_filter() \ 1636Is a version of make_fir_filter; \ 1637it returns one of the standard smoothing filters from \ 1638the era when computers were human beings.") 1639 def make_spencer_filter 1640 data = vct(-3, -6, -5, 3, 21, 46, 67, 74, 67, 46, 21, 3, -5, -6, -3) 1641 data.map! do |n| 1642 n / 320.0 1643 end 1644 make_fir_filter(data.length, data) 1645 end 1646 1647 # any-random 1648 # 1649 # arbitrary random number distributions via the "rejection method" 1650 def any_random(amount, e = false) 1651 if amount.zero? 1652 0.0 1653 else 1654 unless e 1655 random(amount) 1656 else 1657 next_random = lambda do | | 1658 len = e.length 1659 x = random(e[len - 2].to_f) 1660 y = random(1.0) 1661 if y <= envelope_interp(x, e) 1662 x 1663 else 1664 next_random.call 1665 end 1666 end.call 1667 end 1668 end 1669 end 1670 1671 def gaussian_distribution(s) 1672 e = [] 1673 den = 2.0 * s * s 1674 x = 0.0 1675 y = -4.0 1676 21.times do |i| 1677 e.push(exp(-((y * y) / den))) 1678 e.push(x) 1679 x += 0.05 1680 y += 0.4 1681 end 1682 e 1683 end 1684 1685 def pareto_distribution(a) 1686 e = [] 1687 scl = 1.0 ** (a + 1.0) / a 1688 x = 0.0 1689 y = 1.0 1690 21.times do |i| 1691 e.push(scl * (a / y ** (a + 1.0))) 1692 e.push(x) 1693 x += 0.05 1694 y += 0.2 1695 end 1696 e 1697 end 1698 # uniform distribution 1699 # map_channel(lambda do |y| any_random(1.0, [0, 1, 1, 1])) 1700 # mostly toward 1.0 1701 # map_channel(lambda do |y| any_random(1.0, [0, 0, 0.95, 0.1, 1, 1])) 1702 # let(gaussian-distribution(1.0)) do |g| 1703 # map_channel(lambda do |y| any_random(1.0, g)) 1704 # end 1705 # let(pareto-distribution(1.0)) do |g| 1706 # map_channel(lambda do |y| any_random(1.0, g)) 1707 # end 1708 1709 # this is the inverse integration function used by CLM to turn a 1710 # distribution function into a weighting function 1711 def inverse_integrate(dist, data_size = 512, e_size = 50) 1712 first_sum = sum = dist[1].to_f 1713 x0 = dist[0].to_f 1714 x1 = dist[-2].to_f 1715 xincr = (x1 - x0) / e_size.to_f 1716 x = x0 1717 e = make_array(e_size * 2) 1718 0.step(e_size * 2 - 1, 2) do |i| 1719 e[i + 1] = x 1720 e[i] = sum 1721 sum += envelope_interp(x, dist) 1722 x += xincr 1723 end 1724 incr = (e[-2] - first_sum) / (data_size - 1) 1725 x = first_sum - incr 1726 make_vct!(data_size) do 1727 x += incr 1728 envelope_interp(x, e) 1729 end 1730 end 1731 1732 def gaussian_envelope(s) 1733 den = 2.0 * s * s 1734 x = -1.0 1735 y = -4.0 1736 e = make_array(42) 1737 0.step(41, 2) do |i| 1738 e[i] = x 1739 e[i + 1] = exp(-((y * y) / den)) 1740 x += 0.1 1741 y += 0.4 1742 end 1743 e 1744 end 1745 # make_rand(:envelope, gaussian-envelope(1.0)) 1746 1747 # Julius Smith stuff 1748 # 1749 # these are from "Mathematics of the DFT", W3K Pubs 1750 add_help(:channel_mean, 1751 "channel_mean(snd, chn) \ 1752Returns the average of the samples in the given channel: <f,1>/n") 1753 def channel_mean(snd, chn) 1754 sum = 0.0 1755 n = framples(snd, chn) 1756 scan_channel(lambda do |y| 1757 sum += y 1758 false 1759 end, 0, n, snd, chn) 1760 sum / n 1761 end 1762 1763 add_help(:channel_total_energy, 1764 "channel_total_energy(snd, chn) \ 1765Returns the sum of the squares of all the samples in the given channel: <f,f>") 1766 def channel_total_energy(snd, chn) 1767 sum = 0.0 1768 scan_channel(lambda do |y| 1769 sum = sum + y * y 1770 false 1771 end, 0, framples(snd, chn), snd, chn) 1772 sum 1773 end 1774 1775 add_help(:channel_average_power, 1776 "channel_average_power(snd, chn) \ 1777Returns the average power in the given channel: <f,f>/n") 1778 def channel_average_power(snd, chn) 1779 channel_total_energy(snd, chn) / framples(snd, chn) 1780 end 1781 1782 add_help(:channel_rms, 1783 "channel_rms(snd, chn) \ 1784Returns the RMS value of the samples in the given channel: sqrt(<f,f>/n)") 1785 def channel_rms(snd, chn) 1786 sqrt(channel_average_power(snd, chn)) 1787 end 1788 1789 add_help(:channel_variance, 1790 "channel_variance(snd, chn) \ 1791Returns the sample variance in the given channel: <f,f>-((<f,1>/ n)^2") 1792 def channel_variance(snd, chn) 1793 n = framples(snd, chn).to_f 1794 mu = (n / (n - 1.0)) * channel_mean(snd, chn) 1795 p = channel_total_energy(snd, chn) 1796 p - mu * mu 1797 end 1798 1799 add_help(:channel_norm, 1800 "channel_norm(snd, chn) \ 1801Returns the norm of the samples in the given channel: sqrt(<f,f>)") 1802 def channel_norm(snd, chn) 1803 sqrt(channel_total_energy(snd, chn)) 1804 end 1805 1806 add_help(:channel_lp, 1807 "channel_lp(p, snd, chn) \ 1808Returns the Lp norm of the samples in the given channel.") 1809 def channel_lp(lp, snd, chn) 1810 sum = 0.0 1811 scan_channel(lambda do |y| 1812 sum = sum + y.abs ** lp 1813 false 1814 end, 0, framples(snd, chn), snd, chn) 1815 sum ** (1.0 / lp) 1816 end 1817 1818 add_help(:channel_lp_inf, 1819 "channel_lp_inf(snd, chn) \ 1820Returns the maxamp in the given channel (the name is \ 1821just math jargon for maxamp).") 1822 def channel_lp_inf(snd, chn) 1823 mx = 0.0 1824 scan_channel(lambda do |y| 1825 mx = [mx, y.abs].max 1826 false 1827 end, 0, framples(snd, chn), snd, chn) 1828 mx 1829 end 1830 1831 add_help(:channel2_inner_product, 1832 "channel2_inner_product(s1, c1, s2, c2) \ 1833Returns the inner-product of the two channels: <f,g>") 1834 def channel2_inner_product(s1, c1, s2, c2) 1835 sum = 0.0 1836 r1 = make_sampler(0, s1, c1) 1837 r2 = make_sampler(0, s2, c2) 1838 framples(s1, c1).times do |i| 1839 sum = sum + r1.call * r2.call 1840 end 1841 sum 1842 end 1843 1844 add_help(:channel2_angle, 1845 "channel2_angle(s1, c1, s2, c2) \ 1846Treats the two channels as vectors, \ 1847returning the ANGLE between them: acos(<f,g>/(sqrt(<f,f>)*sqrt(<g,g>)))") 1848 def channel2_angle(s1, c1, s2, c2) 1849 inprod = channel2_inner_product(s1, c1, s2, c2) 1850 norm1 = channel_norm(s1, c1) 1851 norm2 = channel_norm(s2, c2) 1852 acos(inprod / (norm1 * norm2)) 1853 end if defined? acos 1854 1855 add_help(:channel2_orthogonal?, 1856 "channel2_orthogonal?(s1, c1, s2, c2) \ 1857Returns true if the two channels' inner-product is 0: <f,g>==0") 1858 def channel2_orthogonal?(s1, c1, s2, c2) 1859 channel2_inner_product(s1, c1, s2, c2).zero? 1860 end 1861 1862 add_help(:channel2_coefficient_of_projection, 1863 "channel2_coefficient_of_projection(s1, c1, s2, c2) \ 1864Returns <f,g>/<f,f>") 1865 def channel2_coefficient_of_projection(s1, c1, s2, c2) 1866 channel2_inner_product(s1, c1, s2, c2) / channel_total_energy(s1, c1) 1867 end 1868 1869 # end of JOS stuff 1870 1871 add_help(:channel_distance, 1872 "channel_distance(s1, c1, s2, c2) \ 1873Returns the euclidean distance between the two channels: sqrt(<f-g,f-g>)") 1874 def channel_distance(s1, c1, s2, c2) 1875 r1 = make_sampler(0, s1, c1) 1876 r2 = make_sampler(0, s2, c2) 1877 sum = 0.0 1878 [framples(s1, c1), framples(s2, c2)].min.times do 1879 diff = r1.call - r2.call 1880 sum = sum + diff * diff 1881 end 1882 sqrt(sum) 1883 end 1884 1885 add_help(:periodogram, 1886 "periodogram(n) \ 1887Displays an N point Bartlett periodogram of the samples in the current channel") 1888 def periodogram(n) 1889 len = framples() 1890 average_data = make_vct(n) 1891 rd = make_sampler(0) 1892 n2 = n * 2 1893 rl = make_vct(n2) 1894 im = make_vct(n2) 1895 len.times do 1896 vct_scale!(rl, 0.0) 1897 vct_scale!(im, 0.0) 1898 n.times do |k| 1899 rl[k] = rd.call 1900 end 1901 mus_fft(rl, im) 1902 n.times do |k| 1903 average_data[k] = average_data[k] + rl[k] * rl[k] + im[k] * im[k] 1904 end 1905 end 1906 graph(vct_scale!(average_data, 1.0 / (len.to_f / n).ceil)) 1907 end 1908 1909 # ssb-am friends 1910 1911 add_help(:shift_channel_pitch, 1912 "shift_channel_pitch(freq, order=40, beg=0, dur=false, \ 1913snd=false, chn=false, edpos=false) \ 1914Uses the ssb-am CLM generator to shift the given channel \ 1915in pitch without changing its length. \ 1916The higher ORDER, the better usually.") 1917 def shift_channel_pitch(freq, order = 40, beg = 0, dur = false, 1918 snd = false, chn = false, edpos = false) 1919 gen = make_ssb_am(freq, order) 1920 map_channel(lambda do |y| 1921 ssb_am(gen, y) 1922 end, 1923 beg, dur, snd, chn, edpos, 1924 format("%s(%s, %s, %s, %s", 1925 get_func_name, freq, order, beg, dur)) 1926 end 1927 1928 add_help(:hz_to_2pi, 1929 "hz_to_2pi(freq) \ 1930Is like hz2radians but uses the current sound's srate, not mus_srate.") 1931 def hz_to_2pi(freq) 1932 (TWO_PI * freq) / srate() 1933 end 1934 1935 def ssb_bank(old_freq, new_freq, pairs, 1936 order = 40, bw = 50.0, beg = 0, dur = false, 1937 snd = false, chn = false, edpos = false) 1938 factor = (new_freq - old_freq.to_f) / old_freq 1939 mx = maxamp 1940 ssbs = make_array(pairs) 1941 bands = make_array(pairs) do |i| 1942 aff = (i + 1.0) * old_freq 1943 bwf = bw * (1.0 + (i + 1.0) / (2.0 * pairs)) 1944 ssbs[i] = make_ssb_am((i + 1.0) * factor * old_freq) 1945 make_bandpass(hz_to_2pi(aff - bwf), hz_to_2pi(aff + bwf), order) 1946 end 1947 as_one_edit_rb("%s(%s, %s, %s, %s, %s, %s, %s", 1948 get_func_name, old_freq, new_freq, 1949 pairs, order, bw, beg, dur) do | | 1950 nmx = 0.0 1951 map_channel_rb(beg, dur, snd, chn, edpos) do |y| 1952 sum = 0.0 1953 ssbs.zip(bands) do |sbs, bds| 1954 sum += ssb_am(sbs, bandpass(bds, y)) 1955 end 1956 nmx = [nmx, sum.abs].max 1957 sum 1958 end 1959 scale_channel(mx / nmx, beg, dur, snd, chn) 1960 end 1961 end 1962 1963 # this version adds a frequency envelope 1964 # ssb_bank_env(557, 880, [0, 0, 1, 100.0], 7) 1965 def ssb_bank_env(old_freq, new_freq, freq_env, pairs, 1966 order = 40, bw = 50.0, beg = 0, dur = false, 1967 snd = false, chn = false, edpos = false) 1968 factor = (new_freq - old_freq.to_f) / old_freq 1969 mx = maxamp 1970 ssbs = make_array(pairs) 1971 frenvs = make_array(pairs) 1972 bands = make_array(pairs) do |i| 1973 aff = (i + 1.0) * old_freq 1974 bwf = bw * (1.0 + (i + 1.0) / (2.0 * pairs)) 1975 ssbs[i] = make_ssb_am((i + 1.0) * factor * old_freq) 1976 frenvs[i] = make_env(:envelope, freq_env, 1977 :scaler, hz2radians(i.to_f), 1978 :length, framples(snd, chn) - 1) 1979 make_bandpass(hz_to_2pi(aff - bwf), hz_to_2pi(aff + bwf), order) 1980 end 1981 as_one_edit_rb("%s(%s, %s, %s, %s, %s, %s, %s, %s", 1982 get_func_name, old_freq, new_freq, freq_env.inspect, 1983 pairs, order, bw, beg, dur) do | | 1984 nmx = 0.0 1985 map_channel_rb(beg, dur, snd, chn, edpos) do |y| 1986 sum = 0.0 1987 ssbs.each_with_index do |sbs, i| 1988 sum += ssb_am(sbs, bandpass(bands[i], y), env(frenvs[i])) 1989 end 1990 nmx = [nmx, sum.abs].max 1991 sum 1992 end 1993 scale_channel(mx / nmx, beg, dur, snd, chn) 1994 end 1995 end 1996 1997 #vct|channel|spectral-polynomial 1998 def vct_polynomial(v, coeffs) 1999 new_v = Vct.new(v.length, coeffs.last) 2000 (coeffs.length - 2).downto(0) do |i| 2001 new_v.multiply!(v).offset!(coeffs[i]) 2002 end 2003 new_v 2004 end 2005 2006 def channel_polynomial(coeffs, snd = false, chn = false) 2007 len = framples(snd, chn) 2008 v = channel2vct(0, len, snd, chn) 2009 vct2channel(vct_polynomial(v, coeffs), 0, len, snd, chn, false, 2010 format("%s(%s", get_func_name, coeffs.to_str)) 2011 end 2012 2013 # channel_polynomial(vct(0.0, 0.5)) = x*0.5 2014 # channel_polynomial(vct(0.0, 1.0, 1.0, 1.0)) = x*x*x + x*x + x 2015 # 2016 # convolution -> * in freq 2017 2018 def spectral_polynomial(coeffs, snd = false, chn = false) 2019 len = framples(snd, chn) 2020 sound = channel2vct(0, len, snd, chn) 2021 num_coeffs = coeffs.length 2022 fft_len = if num_coeffs < 2 2023 len 2024 else 2025 (2.0 ** (log((num_coeffs - 1.0) * len) / log(2.0)).ceil).to_i 2026 end 2027 rl1 = make_vct(fft_len) 2028 rl2 = make_vct(fft_len) 2029 new_sound = make_vct(fft_len) 2030 if coeffs[0] > 0.0 2031 new_sound.map! do 2032 mus_random(coeffs[0]) 2033 end 2034 end 2035 if num_coeffs > 1 2036 new_sound.add!(sound.scale(coeffs[1])) 2037 if num_coeffs > 2 2038 peak = maxamp(snd, chn) 2039 rl1.scale!(0.0).add!(sound) 2040 (2...num_coeffs).each do |i| 2041 convolution(rl1, rl2.scale!(0.0).add(sound), fft_len) 2042 new_sound.add!(rl1.scale((coeffs[i] * peak) / rl1.peak)) 2043 end 2044 new_sound.scale!(peak / new_sound.peak) 2045 end 2046 end 2047 vct2channel(new_sound, 0, [len, len * (num_coeffs - 1)].max, 2048 snd, chn, false, 2049 format("%s(%s", get_func_name, coeffs.to_str)) 2050 end 2051 2052 # SCENTROID 2053 # 2054 # by Bret Battey 2055 # Version 1.0 July 13, 2002 2056 # translated to Snd/Scheme Bill S 19-Jan-05 2057 # 2058 # Returns the continuous spectral centroid envelope of a sound. The 2059 # spectral centroid is the "center of gravity" of the spectrum, and it 2060 # has a rough correlation to our sense of "brightness" of a sound. 2061 # 2062 # [Beauchamp, J., "Synthesis by spectral amplitude and 'brightness' matching 2063 # analyzed musical sounds". 2064 # Journal of Audio Engineering Society 30(6), 396-406] 2065 # 2066 # The formula used is: 2067 # C = [SUM<n=1toj>F(n)A(n)] / [SUM<n=1toj>A(n)] 2068 # Where j is the number of bins in the analysis, 2069 # F(n) is the frequency of a given bin, 2070 # A(n) is the magnitude of the given bin. 2071 # 2072 # If a pitch envelope for the analyzed sound is available, the results 2073 # of SCENTROID can be used with the function NORMALIZE-CENTROID, 2074 # below, to provide a "normalized spectral centroid". 2075 # 2076 # DB-FLOOR -- Frames below this decibel level (0 dB = max) will be 2077 # discarded and returned with spectral centroid = 0 2078 # 2079 # RFREQ -- Rendering frequency. Number of measurements per second. 2080 # 2081 # FFTSIZE -- FFT window size. Must be a power of 2. 4096 is 2082 # recommended. 2083 2084 add_help(:scentroid, 2085 "scentroid(file, beg=0, dur=false, db_floor=-40, \ 2086rfreq=100, fftsize=4096) \ 2087Returns the spectral centroid envelope of a sound; \ 2088RFREQ is the rendering frequency, the number of measurements per second; \ 2089DB_FLOOR is the level below which data will be ignored.") 2090 def scentroid(file, beg = 0.0, dur = false, db_floor = -40, 2091 rfreq = 100.0, fftsize = 4096) 2092 assert_type(File.exist?(file), file, 0, "an existing file") 2093 fsr = srate(file) 2094 incrsamps = (fsr / rfreq).floor 2095 start = (beg * fsr).floor 2096 ende = (start + (dur ? (dur * fsr).floor : (framples(file) - beg))).floor 2097 fdr = make_vct(fftsize) 2098 fdi = make_vct(fftsize) 2099 windows = ((ende - start) / incrsamps).floor + 1 2100 results = make_vct(windows) 2101 fft2 = (fftsize / 2.0).floor 2102 binwidth = fsr / fftsize.to_f 2103 rd = make_readin(file) 2104 loc = 0 2105 start.step(ende, incrsamps) do |i| 2106 rd.location = i 2107 sum_of_squares = 0.0 2108 fdr.map! do 2109 val = readin(rd) 2110 sum_of_squares = sum_of_squares + val * val 2111 val 2112 end 2113 if linear2db(sqrt(sum_of_squares / fftsize.to_f)) >= db_floor 2114 numsum = 0.0 2115 densum = 0.0 2116 fdi.map! do |x| 2117 0.0 2118 end 2119 mus_fft(fdr, fdi, fftsize) 2120 rectangular2polar(fdr, fdi) 2121 fft2.times do |j| 2122 numsum = numsum + j * binwidth * fdr[j] 2123 densum += fdr[j] 2124 end 2125 results[loc] = numsum / densum 2126 end 2127 loc += 1 2128 end 2129 results 2130 end 2131 2132 # invert_filter inverts an FIR filter 2133 # 2134 # say we previously filtered a sound via filter_channel(vct(0.5, 0.25, 0.125)) 2135 # and we want to undo it without using undo_edit: 2136 # filter_channel(invert_filter(vct(0.5, 0.25, 0.125))) 2137 # 2138 # there are a million gotchas here. The primary one is that the inverse 2139 # filtercan "explode" -- the coefficients can grow without bound. For 2140 # example, any filter returned by spectrum2coeffs above will be a problem 2141 # (it always returns a "linear phase" filter). 2142 2143 add_help(:invert_filter, 2144 "invert_filter(coeffs) \ 2145Tries to return an inverse filter to undo the effect of the FIR filter coeffs.") 2146 def invert_filter(fcoeffs) 2147 order = fcoeffs.length + 32 2148 coeffs = Vct.new(order) 2149 fcoeffs.each_with_index do |val, i| 2150 coeffs[i] = val 2151 end 2152 nfilt = Vct.new(order) 2153 nfilt[0] = 1.0 / coeffs.first 2154 (1...order).each do |i| 2155 sum = 0.0 2156 k = i 2157 i.times do |j| 2158 sum = sum + nfilt[j] * coeffs[k] 2159 k -= 1 2160 nfilt[i] = sum / -coeffs.first 2161 end 2162 end 2163 nfilt 2164 end 2165 2166 # Volterra filter 2167 # 2168 # one of the standard non-linear filters 2169 # this version is taken from Monson Hayes "Statistical DSP and Modeling" 2170 # it is a slight specialization of the form mentioned by J O Smith and others 2171 2172 class Volterra_filter < Musgen 2173 def initialize(acoeffs, bcoeffs) 2174 super() 2175 @as = acoeffs 2176 @bs = bcoeffs 2177 @xs = Vct.new([acoeffs.length, bcoeffs.length].max) 2178 end 2179 2180 def inspect 2181 format("%s.new(%s, %s)", self.class, @as.to_str, @bs.to_str) 2182 end 2183 2184 def to_s 2185 format("#<%s acoeffs: %s, bcoeffs: %s>", self.class, @as, @bs) 2186 end 2187 2188 def run_func(val1 = 0.0, val2 = 0.0) 2189 volterra_filter(val1) 2190 end 2191 2192 def volterra_filter(x) 2193 xlen = @xs.length 2194 @xs.move!(xlen - 1, xlen - 2, true) 2195 @xs.first = x 2196 sum = dot_product(@as, @xs, @as.length) 2197 @bs.length.times do |i| 2198 @bs.length.times do |j| 2199 sum = sum + @bs[j] * @xs[i] * @xs[j] 2200 end 2201 end 2202 sum 2203 end 2204 end 2205 2206 add_help(:make_volterra_filter, 2207 "make_volterra_filter(acoeffs, bcoeffs) \ 2208Returns a list for use with volterra-filter, \ 2209producing one of the standard non-linear filters.") 2210 def make_volterra_filter(acoeffs, bcoeffs) 2211 Volterra_filter.new(acoeffs, bcoeffs) 2212 end 2213 2214 add_help(:volterra_filter, 2215 "volterra_filter(flt, x) \ 2216Takes FLT, a Volterra_filter object returned by make_volterra_filter, \ 2217and an input X, and returns the (non-linear filtered) result.") 2218 def volterra_filter(flt, x) 2219 flt.volterra_filter(x) 2220 end 2221 # flt = make_volterra_filter(vct(0.5, 0.1), vct(0.3, 0.2, 0.1)) 2222 # map_channel(lambda do |y| volterra_filter(flt, y) end) 2223 2224 # moving-sum generator (the sum norm or 1-norm) 2225 2226 add_help(:make_moving_sum, 2227 "make_moving_sum(size=128) \ 2228Returns a moving-sum generator. \ 2229The generator keeps a running window of the last SIZE inputs, \ 2230returning the sum of the absolute values of the samples in that window.") 2231 def make_moving_sum(size = 128) 2232 gen = make_moving_average(size) 2233 gen.increment = 1.0 2234 gen 2235 end 2236 2237 add_help(:moving_sum, 2238 "moving_sum(gen, y) \ 2239Returns the sum of the absolute values in a moving \ 2240window over the last few inputs.") 2241 def moving_sum(gen, y) 2242 moving_average(gen, y.abs) 2243 end 2244 2245 def make_unmoving_sum() 2246 make_one_pole(1.0, -1.0) 2247 end 2248 alias unmoving_sum one_pole 2249 2250 # moving-rms generator 2251 2252 add_help(:make_moving_rms, 2253 "make_moving_rms(size=128) \ 2254Returns a moving-rms generator. \ 2255The generator keeps a running window of the last SIZE inputs, \ 2256returning the rms of the samples in that window.") 2257 def make_moving_rms(size = 128) 2258 make_moving_average(size) 2259 end 2260 2261 add_help(:moving_rms, 2262 "moving_rms(gen, y) \ 2263Returns the rms of the values in a window over the last few inputs.") 2264 def moving_rms(gen, y) 2265 sqrt([0.0, moving_average(gen, y * y)].max) 2266 end 2267 2268 # moving-length generator (euclidean norm or 2-norm) 2269 2270 add_help(:make_moving_length, 2271 "make_moving_length(size=128) \ 2272Returns a moving-length generator. \ 2273The generator keeps a running window of the last SIZE inputs, \ 2274returning the euclidean length of the vector in that window.") 2275 alias make_moving_length make_moving_sum 2276 2277 add_help(:moving_length, 2278 "moving_length(gen, y) \ 2279Returns the length of the values in a window over the last few inputs.") 2280 alias moving_length moving_rms 2281 2282 # harmonicizer 2283 # (each harmonic is split into a set of harmonics via Chebyshev polynomials) 2284 # obviously very similar to ssb_bank above, but splits harmonics 2285 # individually, rather than pitch-shifting them 2286 add_help(:harmonicizer, 2287 "harmonicizer(freq, coeffs, pairs, order=40, bw=50.0, \ 2288beg=0, dur=false, snd=false, chn=false, edpos=false) \ 2289Splits out each harmonic and replaces it with the spectrum given in coeffs.") 2290 def harmonicizer(freq, coeffs, pairs, 2291 order = 40, 2292 bw = 50.0, 2293 beg = 0, 2294 dur = false, 2295 snd = false, 2296 chn = false, 2297 edpos = false) 2298 bands = make_array(pairs) 2299 pcoeffs = partials2polynomial(coeffs) 2300 avgs = make_array(pairs) 2301 peaks = make_array(pairs) 2302 flt = make_filter(2, vct(1, -1), vct(0, -0.9)) 2303 old_mx = maxamp 2304 new_mx = 0.0 2305 ctr = 40 2306 1.upto(pairs) do |i| 2307 aff = i * freq 2308 bwf = bw * (1.0 + i / (2 * pairs)) 2309 peaks[i - 1] = make_moving_max(128) 2310 avgs[i - 1] = make_moving_average(128) 2311 bands[i - 1] = make_bandpass(hz_to_2pi(aff - bwf), 2312 hz_to_2pi(aff + bwf), order) 2313 end 2314 as_one_edit_rb do 2315 map_channel_rb(beg, dur, snd, chn, edpos) do |y| 2316 sum = 0.0 2317 bands.zip(peaks, avgs) do |bs, ps, as| 2318 sig = bandpass(bs, y) 2319 mx = moving_max(ps, sig) 2320 amp = moving_average(as, mx > 0.0 ? [100.0, 1.0 / mx].min : 0.0) 2321 if amp > 0.0 2322 sum = sum + mx * polynomial(pcoeffs, amp * sig) 2323 end 2324 end 2325 val = filter(flt, sum) 2326 new_mx = [new_mx, val.abs].max 2327 if ctr.zero? 2328 val 2329 else 2330 ctr -= 1 2331 0.0 2332 end 2333 end 2334 if new_mx > 0.0 2335 scale_channel(old_mx / new_mx, beg, dur, snd, chn) 2336 end 2337 end 2338 end 2339 2340 # linear sampling rate conversion 2341 2342 add_help(:linear_src_channel, 2343 "linear_src_channel(sr, snd=false, chn=false) \ 2344Performs sampling rate conversion using linear interpolation.") 2345 def linear_src_channel(srinc, snd = false, chn = false) 2346 rd = make_sampler(0, snd, chn) 2347 last = rd.call 2348 nxt = rd.call 2349 intrp = 0.0 2350 tempfile = with_sound(:clm, true, 2351 :output, snd_tempnam, 2352 :srate, srate(snd)) do 2353 samp = 0 2354 until sampler_at_end?(rd) 2355 if (pos = intrp) >= 1.0 2356 pos.floor.times do |i| 2357 last, nxt = nxt, rd.call 2358 end 2359 pos -= pos.floor 2360 end 2361 intrp = pos + pos.floor 2362 out_any(samp, last + pos * (nxt - last), 0, $output) 2363 samp += 1 2364 end 2365 end.output 2366 len = mus_sound_framples(tempfile) 2367 set_samples(0, len - 1, tempfile, snd, chn, true, 2368 "%s(%s", get_func_name, srinc, 0, false, true) 2369 # first true=truncate to new length, false=at current edpos, 2370 # true=auto delete temp file 2371 end 2372 2373 # Mathews/Smith High-Q filter as described in 2374 # http://ccrma.stanford.edu/~jos/smac03maxjos/ 2375 2376 class Mfilter < Musgen 2377 def initialize(decay, freq) 2378 super() 2379 @decay = decay.to_f 2380 @frequency = freq.to_f 2381 @eps = 2.0 * sin((PI * freq) / mus_srate()) 2382 @xn = @yn = 0.0 2383 end 2384 attr_accessor :decay, :eps 2385 2386 def inspect 2387 format("%s.new(%0.3f, %0.3f)", self.class, @decay, @frequency) 2388 end 2389 2390 def to_s 2391 format("#<%s decay: %0.3f, frequency: %0.3f>", 2392 self.class, @decay, @frequency) 2393 end 2394 2395 def mfilter(x_input = 0.0, y_input = 0.0) 2396 @xn = x_input + @decay * (@xn - @eps * @yn) 2397 @yn = y_input + @decay * (@eps * @xn + @yn) 2398 end 2399 end 2400 2401 def make_mfilter(*args) 2402 Mfilter.new(get_args(args, :decay, 0.99), 2403 get_args(args, :frequency, 1000.0)) 2404 end 2405 2406 def mfilter(m, x_input = 0.0, y_input = 0.0) 2407 m.mfilter(x_input, y_input) 2408 end 2409=begin 2410 with_sound(:clm, true) do 2411 rd = make_sampler(0, "now.snd") 2412 m = make_mfilter 2413 10000.times do |i| outa(i, mfilter(m, 0.1 * rd.call), $output) end 2414 end 2415=end 2416 # 2417 # sweep center freq: 2418=begin 2419 with_sound(:clm, true) do 2420 rd = make_sampler(0, "oboe.snd") 2421 m = make_mfilter(:decay, 0.99, :frequency, 1000) 2422 e = make_env([0, 100, 1, 2000], :length, 10000) 2423 10000.times do |i| 2424 outa(i, mfilter(m, 0.1 * rd.call), $output) 2425 m.eps = 2.0 * sind((PI * env(e)) / mus_srate()) 2426 end 2427 end 2428=end 2429 # 2430 # harmonics: 2431=begin 2432 with_sound(:clm, true, :statistics, true) do 2433 noi = make_rand(10000) 2434 filters = make_array(9) do 2435 make_mfilter(:decay, 0.999, :frequency, 400.0 * (i + 1)) 2436 end 2437 10000.times do |i| 2438 sum = 0.0 2439 input = 0.01 * rand(noi) 2440 filters.each do |f| sum = sum + (1.0 / (j + 1)) * mfilter(f, input) end 2441 outa(i, sum $output) 2442 end 2443 end 2444=end 2445 2446 # 2447 # spectrum displayed in various frequency scales 2448 # 2449 class Display_bark_fft 2450 # click in lisp-graph to change the tick placement choice 2451 def initialize 2452 @bark_fft_size = 0 2453 @bark_tick_function = 0 2454 end 2455 attr_reader :bark_tick_function 2456 2457 def display_bark_fft(snd, chn) 2458 ls = left_sample(snd, chn) 2459 rs = right_sample(snd, chn) 2460 fftlen = (2 ** (log((rs - ls) + 1.0) / log(2.0)).ceil.to_i).to_i 2461 if fftlen > 0 2462 data = channel2vct(ls, fftlen, snd, chn) 2463 normalized = (transform_normalization(snd, chn) != Dont_normalize) 2464 linear = true 2465 if vct?(data) 2466 fft = snd_spectrum(data, 2467 fft_window(snd, chn), fftlen, linear, 2468 fft_window_beta(snd, chn), false, normalized) 2469 if vct?(fft) 2470 sr = srate(snd) 2471 mx = fft.peak 2472 data_len = fft.length 2473 2474 # bark settings 2475 bark_low = bark(20.0).floor 2476 bark_high = bark(0.5 * sr).ceil 2477 bark_frqscl = data_len / (bark_high - bark_low) 2478 bark_data = Vct.new(data_len) 2479 2480 # mel settings 2481 mel_low = mel(20.0).floor 2482 mel_high = mel(0.5 * sr).ceil 2483 mel_frqscl = data_len / (mel_high - mel_low) 2484 mel_data = Vct.new(data_len) 2485 2486 # erb settings 2487 erb_low = erb(20.0).floor 2488 erb_high = erb(0.5 * sr).ceil 2489 erb_frqscl = data_len / (erb_high - erb_low) 2490 erb_data = Vct.new(data_len) 2491 2492 @bark_fft_size = fftlen 2493 fftlenf = fftlen.to_f 2494 2495 fft.each_with_index do |val, i| 2496 frq = sr * (i / fftlenf) 2497 bark_bin = (bark_frqscl * (bark(frq) - bark_low)).round 2498 mel_bin = (mel_frqscl * (mel(frq) - mel_low)).round 2499 erb_bin = (erb_frqscl * (erb(frq) - erb_low)).round 2500 if bark_bin.between?(0, data_len - 1) 2501 bark_data[bark_bin] += val 2502 end 2503 if mel_bin.between?(0, data_len - 1) 2504 mel_data[mel_bin] += val 2505 end 2506 if erb_bin.between?(0, data_len - 1) 2507 erb_data[erb_bin] += val 2508 end 2509 end 2510 2511 if normalized 2512 bmx = bark_data.peak 2513 mmx = mel_data.peak 2514 emx = erb_data.peak 2515 if (mx - bmx).abs > 0.01 2516 bark_data.scale!(mx / bmx) 2517 end 2518 if (mx - mmx).abs > 0.01 2519 mel_data.scale!(mx / mmx) 2520 end 2521 if (mx - emx).abs > 0.01 2522 erb_data.scale!(mx / emx) 2523 end 2524 end 2525 graph([bark_data, mel_data, erb_data], 2526 "ignored", 2527 20.0, 0.5 * sr, 2528 0.0, (normalized ? 1.0 : data_len * y_zoom_slider(snd, chn)), 2529 snd, chn, 2530 false, Show_bare_x_axis) 2531 end 2532 end 2533 end 2534 false 2535 end 2536 2537 def mark_bark_labels(snd, chn) 2538 # at this point the x axis has no markings, but there is room 2539 # for labels and ticks 2540 old_foreground_color = foreground_color(snd, chn, Copy_context) 2541 # assume at start the foreground color is correct 2542 axinfo = axis_info(snd, chn, Lisp_graph) 2543 axis_x0 = axinfo[10] 2544 axis_x1 = axinfo[12] 2545 axis_y0 = axinfo[13] 2546 axis_y1 = axinfo[11] 2547 label_height = 15 2548 char_width = 8 2549 sr2 = 0.5 * srate(snd) 2550 minor_tick_len = 6 2551 major_tick_len = 12 2552 tick_y0 = axis_y1 2553 minor_y0 = axis_y1 + minor_tick_len 2554 major_y0 = axis_y1 + major_tick_len 2555 bark_label_font = snd_font(3) 2556 bark_numbers_font = snd_font(2) 2557 label_pos = (axis_x0 + 0.45 * (axis_x1 - axis_x0)).to_i 2558 cr = channel_widgets(snd, chn)[17] 2559 scale_position = lambda do |scale, f| 2560 b20 = scale.call(20.0) 2561 (axis_x0 + 2562 ((axis_x1 - axis_x0) * (scale.call(f) - b20)) / 2563 (scale.call(sr2) - b20)).round 2564 end 2565 bark_position = lambda do |f| 2566 scale_position.call(method(:bark).to_proc, f) 2567 end 2568 mel_position = lambda do |f| 2569 scale_position.call(method(:mel).to_proc, f) 2570 end 2571 erb_position = lambda do |f| 2572 scale_position.call(method(:erb).to_proc, f) 2573 end 2574 draw_bark_ticks = lambda do |bark_function| 2575 if bark_numbers_font 2576 set_current_font(bark_numbers_font, snd, chn, Copy_context) 2577 end 2578 draw_line(axis_x0, tick_y0, axis_x0, major_y0, 2579 snd, chn, Copy_context, cr) 2580 i1000 = scale_position.call(bark_function, 1000.0) 2581 i10000 = scale_position.call(bark_function, 10000.0) 2582 draw_line(i1000, tick_y0, i1000, major_y0, 2583 snd, chn, Copy_context, cr) 2584 draw_line(i10000, tick_y0, i10000, major_y0, 2585 snd, chn, Copy_context, cr) 2586 draw_string("20", axis_x0, major_y0, 2587 snd, chn, Copy_context, cr) 2588 draw_string("1000", i1000 - 3 * 4, major_y0, 2589 snd, chn, Copy_context, cr) 2590 draw_string("10000", i10000 - 6 * 4, major_y0, 2591 snd, chn, Copy_context, cr) 2592 draw_string("fft size: #{@bark_fft_size}", axis_x0 + 10, axis_y0, 2593 snd, chn, Copy_context, cr) 2594 100.step(1000, 100) do |i| 2595 i100 = scale_position.call(bark_function, i) 2596 draw_line(i100, tick_y0, i100, minor_y0, snd, chn, Copy_context, cr) 2597 end 2598 2000.step(10000, 1000) do |i| 2599 i1000 = scale_position.call(bark_function, i) 2600 draw_line(i1000, tick_y0, i1000, minor_y0, snd, chn, Copy_context, cr) 2601 end 2602 end 2603 2604 # bark label/ticks 2605 if self.bark_tick_function.zero? 2606 draw_bark_ticks.call(bark_position) 2607 end 2608 if bark_label_font 2609 set_current_font(bark_label_font, snd, chn, Copy_context) 2610 end 2611 draw_string("bark,", label_pos, axis_y1 + label_height, 2612 snd, chn, Copy_context, cr) 2613 # mel label/ticks 2614 set_foreground_color(snd_color(2), snd, chn, Copy_context) 2615 if self.bark_tick_function == 1 2616 draw_bark_ticks.call(mel_position) 2617 end 2618 if bark_label_font 2619 set_current_font(bark_label_font, snd, chn, Copy_context) 2620 end 2621 draw_string("mel,", char_width * 6 + label_pos, axis_y1 + label_height, 2622 snd, chn, Copy_context, cr) 2623 # erb label/ticks 2624 set_foreground_color(snd_color(4), snd, chn, Copy_context) 2625 if self.bark_tick_function == 2 2626 draw_bark_ticks.call(erb_position) 2627 end 2628 if bark_label_font 2629 set_current_font(bark_label_font, snd, chn, Copy_context) 2630 end 2631 draw_string("erb", 2632 char_width * (6 + 5) + label_pos, 2633 axis_y1 + label_height, 2634 snd, chn, Copy_context, cr) 2635 set_foreground_color(old_foreground_color, snd, chn, Copy_context) 2636 end 2637 2638 # mouse click = move to next scale's ticks 2639 def choose_bark_ticks(snd, chn, button, state, x, y, axis) 2640 if axis == Lisp_graph 2641 @bark_tick_function += 1 2642 if @bark_tick_function > 2 2643 @bark_tick_function = 0 2644 end 2645 update_lisp_graph(snd, chn) 2646 end 2647 end 2648 2649 private 2650 def bark(f) 2651 f2 = f / 7500.0 2652 13.5 * atan(0.00076 * f) + (3.5 * atan(f2 * f2)) 2653 end 2654 2655 def mel(f) 2656 1127.0 * log(1.0 + f / 700.0) 2657 end 2658 2659 def erb(f) 2660 43.0 + 11.17 * log((f + 312.0) / (f + 14675.0)) 2661 end 2662 end 2663 2664 # user's view of display-bark-fft function 2665 def display_bark_fft(off = false) 2666 if off.kind_of?(FalseClass) 2667 db = Display_bark_fft.new 2668 $lisp_graph_hook.add_hook!("display-bark-fft") do |snd, chn| 2669 db.display_bark_fft(snd, chn) 2670 end 2671 $after_lisp_graph_hook.add_hook!("make-bark-label") do |snd, chn| 2672 db.mark_bark_labels(snd, chn) 2673 end 2674 $mouse_click_hook.add_hook!("choose-bark-ticks") do |s, c, b, st, x, y, a| 2675 db.choose_bark_ticks(s, c, b, st, x, y, a) 2676 end 2677 Snd.sounds.each do |snd| 2678 channels(snd).times do |chn| 2679 update_lisp_graph(snd, chn) 2680 end 2681 end 2682 else 2683 $lisp_graph_hook.remove_hook!("display-bark-fft") 2684 $after_lisp_graph_hook.remove_hook!("make-bark-label") 2685 $mouse_click_hook.remove_hook!("choose-bark-ticks") 2686 Snd.sounds.each do |snd| 2687 channels(snd).times do |chn| 2688 set_lisp_graph?(false, snd, chn) 2689 end 2690 end 2691 end 2692 off 2693 end 2694 2695 def undisplay_bark_fft 2696 display_bark_fft(true) 2697 end 2698end 2699 2700include Dsp 2701 2702# dsp.rb ends here 2703