# dsp.rb -- dsp.scm --> dsp.rb # Translator: Michael Scholz # Created: 2005/03/07 13:50:44 # Changed: 2018/04/15 22:52:00 # comments are taken mostly from dsp.scm # # module Dsp # src_duration(e) # dolph(n, gamma) # dolph_1(n, gamma) # down_oct(n, snd, chn) # edot_product(freq, data) # stretch_sound_via_dft(factor, snd, chn) # compute_uniform_circular_string(size, x0, x1, x2, mass, xspring, damp) # testunif(mass, xspring, damp) # test_scanned_synthesis(amp, dur, mass, xspring, damp) # compute_string(size, x0, x1, x2, masses, xsprings, esprings, damps, haptics) # freqdiv(n, snd, chn) # adsat(size, beg, dur, snd, chn) # spike(snd, chn) # spot_freq(samp, snd, chn) # class Flanger < Musgen # initialize(time, amount, speed) # inspect # to_s # run_func(val1, val2) # flanger(inval) # make_flanger(time, amount, speed) # flanger?(obj) # flanger(gen, inval) # chorus(size) # chordalize(amount, base, chord) # zero_phase(snd, chn) # rotate_phase(func, snd, chn) # class Asyfm < Musgen # initialize(*args) # inspect # to_s # run_func(val1, val2) # asyfm_J(input) # asyfm_I(input) # make_asyfm(*args) # asyfm?(obj) # asyfm_J(gen, input) # asyfm_I(gen, input) # cosine_summation(gen, r) # kosine_summation(gen, r, k) # fejer_sum(angle, n) # legendre_sum(angle, n) # sum_of_n_sines(angle, n) # sum_of_n_odd_sines(angle, n) # sum_of_n_odd_cosines(angle, n) # band_limited_sawtooth(x, a, n, fi) # brighten_slightly(amount, snd, chn) # brighten_slightly_1(coeffs, snd, chn) # spectrum2coeffs(order, spectr) # fltit_1(order, spectr) # make_hilbert_transform(len) # make_highpass(fc, len) # make_lowpass(fc, len) # make_bandpass(flo, fhi, len) # make_bandstop(flo, fhi, len) # make_differentiator(len) # make_butter_high_pass(freq) # make_butter_low_pass(freq) # make_butter_band_pass(freq, band) # make_butter_band_reject(freq, band) # make_biquad(a0, a1, a2, b1, b2) # make_iir_low_pass_2(fc, din) # make_iir_high_pass_2(fc, din) # make_iir_band_pass_2(f1, f2) # make_iir_band_stop_2(f1, f2) # make_eliminate_hum(hum_freq, hum_harmonics, bandwidth) # eliminate_hum(gens, x0) # make_peaking_2(f1, f2, m) # cascade2canonical(a) # make_butter_lp(m, fc) # make_butter_hp(m, fc) # make_butter_bp(m, f1, f2) # make_butter_bs(m, f1, f2) # make_notch_frequency_response(cur_srate, freqs, notch_width) # notch_channel(freqs, filter_order, beg, dur, snd, chn, edpos, trc, n_width) # notch_sound(freqs, filter_order, snd, chn, notch_width) # notch_selection(freqs, filter_order, snd, chn, notch_width) # fractional_fourier_transform(fr, fi, n, v) # z_transform(f, n, z) # dht(data) # find_sine(freq, beg, dur, snd) # goertzel(freq, beg, dur, snd) # make_spencer_filter # any_random(amount, e) # gaussian_distribution(s) # pareto_distribution(a) # inverse_integrate(dist, data_size, e_size) # gaussian_envelope(s) # channel_mean(snd, chn) # channel_total_energy(snd, chn) # channel_average_power(snd, chn) # channel_rms(snd, chn) # channel_variance(snd, chn) # channel_norm(snd, chn) # channel_lp(p, snd, chn) # channel_lp_inf(snd, chn) # channel2_inner_product(s1, c1, s2, c2) # channel2_angle(s1, c1, s2, c2) # channel2_orthogonal?(s1, c1, s2, c2) # channel2_coefficient_of_projection(s1, c1, s2, c2) # channel_distance(s1, c1, s2, c2) # periodogram(n) # shift_channel_pitch(freq, order, beg, dur, snd, chn, edpos) # hz_to_2pi(freq) # ssb_bank(old_freq, new_freq, pairs, order, bw, beg, dur, snd, chn, edpos) # ssb_bank_env(old_freq, new_freq, fq_env, pairs, ord, bw, beg, dur, s, c, ep) # # vct_polynomial(v, coeffs) # channel_polynomial(coeffs, snd, chn) # spectral_polynomial(coeffs, snd, chn) # scentroid(file, beg, dur, db_floor, rfreq, fftsize) # invert_filter(fcoeffs) # # class Volterra_filter < Musgen # initialize(acoeffs, bcoeffs) # inspect # to_s # run_func(val1, val2) # volterra_filter(x) # # make_volterra_filter(acoeffs, bcoeffs) # volterra_filter(flt, x) # make_moving_sum(size) # moving_sum(gen, y) # make_moving_rms(size) # moving_rms(gen, y) # make_moving_length(size) # moving_length(gen, y) # harmonicizer(freq, coeffs, pairs, order, bw, beg, dur, snd, chn, edpos) # linear_src_channel(srinc, snd, chn) # # class Mfilter < Musgen # initialize(decay, freq) # inspect # to_s # mfilter(x_input, y_input) # # make_mfilter(*args) # mfilter(m, x_input, y_input) # # class Display_bark_fft # display_bark_fft(snd, chn) # mark_bark_labels(snd, chn) # choose_bark_ticks(snd, chn, button, state, x, y, axis) # # display_bark_fft(off) # undisplay_bark_fft require "ws" require "env" module Dsp # src_duration (see src-channel in extsnd.html) add_help(:src_duration, "src_duration(envelope) \ Returns the new duration of a sound after using ENVELOPE for \ time-varying sampling-rate conversion.") def src_duration(e) e.map! do |x| x.to_f end ex0 = e.first ex1 = e[-2] all_x = ex1 - ex0 dur = 0.0 0.step(e.length - 3, 2) do |i| x0, xy0, x1, xy1 = e[i, 4] y0 = xy0.zero? ? 1.0 : (1.0 / xy0) y1 = xy1.zero? ? 1.0 : (1.0 / xy1) area = if (xy0 - xy1).abs < 0.0001 y0 * ((x1 - x0) / all_x) else ((log(y1) - log(y0)) / (xy0 - xy1)) * ((x1 - x0) / all_x) end dur += area.abs end dur end # Dolph-Chebyshev window # # formula taken from Richard Lyons, "Understanding DSP" # see clm.c for C version (using either GSL's or GCC's complex trig functions) add_help(:dolph, "dolph(n, gamma) \ Produces a Dolph-Chebyshev FFT data window of N points \ using GAMMA as the window parameter.") def dolph(n, gamma) alpha = cosh(acosh(10.0 ** gamma) / n) den = 1.0 / cosh(n * acosh(alpha)) freq = PI / n rl = make_vct(n) im = make_vct(n) phase = 0.0 n.times do |i| val = den * cos(n * acos(alpha * cos(phase))) rl[i] = val.real im[i] = val.imag phase += freq end fft(rl, im, -1) vct_scale!(rl, 1.0 / vct_peak(rl)) j = n / 2 n.times do |i| im[i] = rl[j] j += 1 if j == n j = 0 end end im end if defined? acosh # this version taken from Julius Smith's "Spectral Audio..." with # three changes it does the DFT by hand, and is independent of # anything from Snd (fft, vcts etc) add_help(:dolph_1, "dolph_1(n, gamma) \ Produces a Dolph-Chebyshev FFT data window of N points \ using GAMMA as the window parameter.") def dolph_1(n, gamma) alpha = cosh(acosh(10.0 ** gamma) / n) den = 1.0 / cosh(n * acosh(alpha)) freq = PI / n vals = make_array(n) w = make_array(n) pk = 0.0 mult = -1.0 phase = -HALF_PI n.times do |i| vals[i] = mult * den * cos(n * acos(alpha * cos(phase))) mult *= -1.0 phase += freq end n.times do |i| sum = 0.0 n.times do |j| sum = sum + vals[j] * exp((2.0 * Complex(0.0, 1.0) * PI * j * i) / n) end w[i] = sum.abs if w[i] > pk pk = w[i] end end w.map! do |val| val / pk end end if defined? acosh # move sound down by n (a power of 2) # I think this is "stretch" in DSP jargon -- to interpolate in the # time domain we're squeezing the frequency domain the power-of-2 # limitation is based on the underlying fft function's insistence on # power-of-2 data sizes see stretch-sound-via-dft below for a # general version add_help(:down_oct, "down_oct(n, snd=false, chn=false) \ Moves a sound down by power of 2 N.") def down_oct(n, snd = false, chn = false) len = framples(snd, chn) pow2 = (log(len) / log(2)).ceil fftlen = (2 ** pow2).round fftscale = 1.0 / fftlen rl1 = channel2vct(0, fftlen, snd, chn) im1 = make_vct(fftlen) fft(rl1, im1, 1) vct_scale!(rl1, fftscale) vct_scale!(im1, fftscale) rl2 = make_vct(2 * fftlen) im2 = make_vct(2 * fftlen) k = fftlen - 1 j = fftlen * n - 1 (0...(fftlen / 2)).each do |i| vct_set!(rl2, i, rl1[i]) vct_set!(rl2, j, rl1[k]) vct_set!(im2, i, im1[i]) vct_set!(im2, j, im1[k]) k -= 1 j -= 1 end fft(rl2, im2, -1) vct2channel(rl2, 0, n * len, snd, chn, false, format("%s(%s", get_func_name, n)) end add_help(:edot_product, "edot_product(freq, data) \ Sum of (e^freq*i) * data[i]") def edot_product(freq, data) sum = 0.0 data.each_with_index do |val, i| sum = sum + exp(i.to_f * freq) * val end sum end unless defined? edot_product add_help(:stretch_sound_via_dft, "stretch_sound_via_dft(factor, snd=false, chn=false) \ Makes the given channel longer (FACTOR should be > 1.0) by \ squeezing in the frequency domain, \ then using the inverse DFT to get the time domain result.") def stretch_sound_via_dft(factor, snd = false, chn = false) factor = factor.to_f n = framples(snd, chn) n2 = (n / 2.0).floor out_n = (n * factor).round in_data = channel2vct(0, n, snd, chn) out_data = make_vct(out_n) fr = make_array(out_n, 0.0) freq = (PI * 2) / n n.times do |i| d = edot_product(freq * Complex(0.0, 1.0) * i, in_data) if i < n2 fr[i] = d else fr[i + (out_n - n - 1)] = d end end freq = (PI * 2) / out_n out_n.times do |i| out_data[i] = (edot_product(freq * Complex(0.0, 1.0) * i, fr) / n).real end vct2channel(out_data, 0, out_n, snd, chn, false, format("%s(%s", get_func_name, factor)) end # compute-uniform-circular-string # # this is a simplification of the underlying table-filling routine # for "scanned synthesis". To watch the wave, open some sound (so # Snd has some place to put the graph), turn off the time domain # display (to give our graph all the window -- to do this in a much # more elegant manner, see snd-motif.scm under scanned-synthesis). def compute_uniform_circular_string(size, x0, x1, x2, mass, xspring, damp) circle_vct_ref = lambda do |v, i| if i < 0 v[i + size] elsif i >= size v[i - size] else v[i] end end dm = damp / mass.to_f km = xspring / mass.to_f denom = 1.0 + dm p1 = (2.0 + (dm - 2.0 * km)) / denom p2 = km / denom p3 = -1.0 / denom size.times do |i| x0[i] = p1 * x1[i] + p2 * (circle_vct_ref.call(x1, i - 1) + circle_vct_ref.call(x1, i + 1)) + p3 * x2[i] end vct_fill!(x2, 0.0) vct_add!(x2, x1) vct_fill!(x1, 0.0) vct_add!(x1, x0) end def testunif(mass, xspring, damp) size = 128 x0 = make_vct(size) x1 = make_vct(size) x2 = make_vct(size) 12.times do |i| x1[i + size / 4 - 6] = sin((TWO_PI * i) / 12.0) end 1024.times do |i| compute_uniform_circular_string(size, x0, x1, x2, mass, xspring, damp) graph(x0, "string", 0, 1.0, -10.0, 10.0) end end def test_scanned_synthesis(amp, dur, mass, xspring, damp) size = 256 x0 = make_vct(size) gx1 = make_vct(size) gx2 = make_vct(size) 12.times do |i| gx1[i + size / 4 - 6] = sin((TWO_PI * i) / 12.0) end gen1 = make_table_lookup(440.0, :wave, gx1) gen2 = make_table_lookup(440.0, :wave, gx2) x1 = gen1.data x2 = gen2.data recompute_samps = 30.0 k = 0.0 kincr = 1.0 / recompute_samps data = make_vct!(dur) do |i| if k >= 1.0 k = 0.0 compute_uniform_circular_string(size, x0, x1, x2, mass, xspring, damp) else k += kincr end g1 = table_lookup(gen1) g2 = table_lookup(gen2) g2 + k * (g1 - g2) end vct_scale!(data, amp / vct_peak(data)) vct2channel(data, 0, dur) end # this is the more general form def compute_string(size, x0, x1, x2, masses, xsprings, esprings, damps, haptics) circle_vct_ref = lambda do |v, i| if i < 0 v[i + size] elsif i >= size v[i - size] else v[i] end end size.times do |i| dm = damps[i] / masses[i] km = xsprings[i] / masses[i] cm = esprings[i] / masses[i] denom = 1.0 + dm + cm p1 = (2.0 + (dm - 2.0 * km)) / denom p2 = km / denom p3 = -1.0 / denom p4 = haptics / (masses[i] * denom) x0[i] = p1 * x1[i] + p2 * (circle_vct_ref.call(x1, i - 1) + circle_vct_ref.call(x1, i + 1)) + p3 * x2[i] + p4 end size.times do |i| x2[i], x1[i] = x1[i], x0[i] end end # "frequency division" -- an effect from sed_sed@my-dejanews.com add_help(:freqdiv, "freqdiv(n, snd=false, chn=false) \ Repeats each nth sample N times (clobbering the intermediate samples): \ freqdiv(8)") def freqdiv(n, snd = false, chn = false) div = 0 curval = 0.0 map_channel(lambda do |val| curval = val if div.zero? div += 1 div = 0 if div == n curval end, 0, false, snd, chn, false, format("%s(%s", get_func_name, n)) end # "adaptive saturation" -- an effect from sed_sed@my-dejanews.com # # a more extreme effect is "saturation": # (map-channel (lambda (val) # (if (< (abs val) .1) val (if (>= val 0.0) 0.25 -0.25)))) add_help(:adsat, "adsat(size, beg=false, dur=false, snd=false, chn=false) \ Is an 'adaptive saturation' sound effect.") def adsat(size, beg = false, dur = false, snd = false, chn = false) mn = 0.0 mx = 0.0 n = 0 vals = make_vct(size) map_channel(lambda do |val| if n == size size.times do |i| if vals[i] >= 0.0 vals[i] = mx else vals[i] = mn end end n = 0 mx = 0.0 mn = 0.0 vals else vals[n] = val mx = val if val > mx mn = val if val < mn n += 1 false end end, beg, dur, snd, chn, false, format("%s(%s, %s, %s", get_func_name, size, beg, dur)) end # spike # # makes sound more spikey -- sometimes a nice effect add_help(:spike, "spike(snd=false, chn=false) \ Multiplies successive samples together to make a sound more spikey.") def spike(snd = false, chn = false) x1 = x2 = 0.0 amp = maxamp(snd, chn) map_channel(lambda do |x0| res = (x0 / (amp * amp)) * x2.abs * x1.abs x2, x1 = x1, x0 res end, 0, false, snd, chn, false, "spike(") end # easily-fooled autocorrelation-based pitch tracker add_help(:spot_freq, "spot_freq(samp=0, snd=false, chn=false) \ Tries to determine the current pitch: spot_freq(left_sample)") def spot_freq(samp = 0, snd = false, chn = false) pow2 = (log(srate(snd) / 20.0) / log(2)).ceil fftlen = (2 ** pow2).round data = autocorrelate(channel2vct(samp, fftlen, snd, chn)) cor_peak = vct_peak(data) cor_peak2 = 2.0 * cor_peak ret = 0.0 (1...fftlen - 2).each do |i| if data[i] < data[i + 1] and data[i + 2] < data[i + 1] logla = log10((cor_peak + data[i]) / cor_peak2) logca = log10((cor_peak + data[i + 1]) / cor_peak2) logra = log10((cor_peak + data[i + 2]) / cor_peak2) offset = (0.5 * (logla - logra)) / (logla + logra + -2.0 * logca) ret = srate(snd) / (2.0 * (i + 1 + offset)) break end end ret end # $graph_hook.add_hook!("examp-left-sample-hook") do |snd, chn, y0, y1| # msg = format("(freq: %.3f)", spot_freq(left_sample(snd, chn))) # status_report(msg, snd) # end # # or # # $mouse_click_hook.add_hook!("examp-cursor-hook") do |snd, chn, # button, state, # x, y, axis| # if axis == Time_graph # status_report(format("(freq: %.3f)", spot_freq(cursor(snd, chn))), snd) # end # end # chorus (doesn't always work and needs speedup) class Flanger < Musgen def initialize(time = 0.05, amount = 20.0, speed = 10.0) super() @time = time @amount = amount @speed = speed @randind = make_rand_interp(:frequency, speed, :amplitude, amount) @data = @randind.data @length = @randind.length len = random(3.0 * time * mus_srate()).floor @flanger = make_delay(:size, len, :max_size, (len + amount + 1).to_i) end def inspect format("%s.new(%s, %s, %s)", self.class, @time, @amount, @speed) end def to_s format("#<%s time: %1.3f, amount: %1.3f, speed: %1.3f>", self.class, @time, @amount, @speed) end def run_func(val1 = 0.0, val2 = 0.0) flanger(val1) end def flanger(inval) inval + delay(@flanger, inval, rand_interp(@randind)) end end def make_flanger(time = 0.05, amount = 20.0, speed = 10.0) Flanger.new(time, amount, speed) end def flanger?(obj) obj.kind_of?(Flanger) end def flanger(gen, inval) gen.flanger(inval) end add_help(:chorus, "chorus(size=5) \ Tries to produce the chorus sound effect.") def chorus(size = 5) dlys = make_array(size) do make_flanger end sum = 0.0 lambda do |inval| dlys.each do |dly| sum += dly.flanger(inval) end sum * 0.25 end end # chordalize (comb filters to make a chord using chordalize-amount # and chordalize-base) add_help(:chordalize, "chordalize(amount=0.95, base=100, chord=[1.00, 0.75, 1.25]) \ Uses harmonically-related comb-filters to bring out a chord in a sound.") def chordalize(amount = 0.95, base = 100, chord = [1.00, 0.75, 1.25]) combs = chord.map do |interval| make_comb(:scaler, amount, :size, (base * interval).round) end scaler = 0.5 / chord.length lambda do |x| val = 0.0 combs.each do |c| val += comb(c, x) end scaler * val end end # zero-phase, rotate-phase # fft games (from the "phazor" package of Scott McNab) add_help(:zero_phase, "zero_phase(snd=false, chn=false) \ Calls fft, sets all phases to 0, and un-ffts.") def zero_phase(snd = false, chn = false) len = framples(snd, chn) pow2 = (log(len) / log(2)).ceil fftlen = (2 ** pow2).round fftscale = 1.0 / fftlen rl = channel2vct(0, fftlen, snd, chn) old_pk = vct_peak(rl) im = make_vct(fftlen) fft(rl, im, 1) rectangular2polar(rl, im) vct_scale!(rl, fftscale) vct_scale!(im, 0.0) fft(rl, im, -1) pk = vct_peak(rl) vct2channel(rl.scale(old_pk / pk), 0, len, snd, chn, false, "zero_phase(") end # (set_)edit_list_proc_counter is defined in clm.rb # It's necessary to produce a uniq method name. add_help(:rotate_phase, "rotate_phase(func, snd=false, chn=false) \ Calls fft, applies func to each phase, then un-ffts.") def rotate_phase(func, snd = false, chn = false) func_name = format("%s_%d", get_func_name, set_edit_list_proc_counter).intern # Proc converted to Method (ie. normal function) for edit_list2function func.to_method(func_name) len = framples(snd, chn) pow2 = (log(len) / log(2)).ceil fftlen = (2 ** pow2).round fftlen2 = (fftlen / 2).floor fftscale = 1.0 / fftlen rl = channel2vct(0, fftlen, snd, chn) im = make_vct(fftlen) old_pk = rl.peak fft(rl, im, 1) rectangular2polar(rl, im) rl.scale!(fftscale) im[0] = 0.0 j = fftlen - 1 (1...fftlen2).each do |i| im[i] = snd_func(func_name, im[i]) im[j] = -im[i] j -= 1 end polar2rectangular(rl, im) fft(rl, im, -1) pk = rl.peak vct2channel(rl.scale(old_pk / pk), 0, len, snd, chn, false, format("%s(Proc.new {|val_r| %s(val_r) }", get_func_name, func_name)) end # rotate_phase(lambda {|x| 0.0 }) # is the same as (zero-phase) # rotate_phase(lambda {|x| random(PI) }) # randomizes phases # rotate_phase(lambda {|x| x }) # returns original # rotate_phase(lambda {|x| -x }) # reverses original # # (might want to write fftlen samps here) # asymmetric FM (bes-i0 case) class Asyfm < Musgen def initialize(*args) super() frequency, ratio, r, index, phase = nil optkey(args, binding, [:frequency, $clm_default_frequency], [:ratio, 1.0], [:r, 1.0], [:index, 1.0], [:phase, 0.0]) @frequency = frequency.to_f @ratio = ratio.to_f @r = r.to_f @index = index.to_f @freq = hz2radians(@frequency) @phase = phase.to_f end attr_accessor :ratio, :r, :index def inspect format("%s.new(:frequency, %s, :ratio, %s, :r, %s, \ :index, %s, :freq, %s, :phase, %s)", self.class, @frequency, @ratio, @r, @index, @freq, @phase) end def to_s format("#<%s freq: %1.3f, ratio: %1.3f, r: %1.3f, \ index: %1.3f, freq: %1.3f, phase: %1.3f>", self.class, @frequency, @ratio, @r, @index, @freq, @phase) end def run_func(val1 = 0.0, val2 = 0.0) asyfm_J(val1) end def asyfm_J(input) # It follows now asyfm-J in generators.scm, not dsp-asyfm-J in clm23.scm. r1 = 1.0 / @r one = ((@r > 1.0) or (@r < 0.0 and @r > -1.0)) ? -1.0 : 1.0 modphase = @ratio * @phase result = exp(0.5 * @index * (@r - r1) * (one + cos(modphase))) * cos(@phase + 0.5 * @index * (@r + r1) * sin(modphase)) @phase = @phase + input + @freq result end def asyfm_I(input) r1 = 1.0 / @r modphase = @ratio * @phase result = exp(0.5 * @index * (@r + r1) * (cos(modphase) - 1.0)) - cos(@phase + 0.5 * @index * (@r - r1) * sin(modphase)) @phase = @phase + input + @freq result end end def make_asyfm(*args) Asyfm.new(*args) end def asyfm?(obj) obj.kind_of?(Asyfm) end def asyfm_J(gen, input) gen.asyfm_J(input) end def asyfm_I(gen, input) gen.asyfm_I(input) end # cosine-summation (a simpler version of sine-summation) # # from Andrews, Askey, Roy "Special Functions" 5.1.16 add_help(:cosine_summation, "cosine_summation(gen, r) \ Is a variant of the CLM sine-summation generator; \ R controls successive sinusoid amplitudes.") def cosine_summation(gen, r) rr = r * r rrp1 = 1.0 + rr rrm1 = 1.0 - rr r2 = 2.0 * r ((rrm1 / (rrp1 - r2 * oscil(gen))) - 1.0) * ((1.0 - r) / r2) end alias make_cosine_summation make_oscil # kosine-summation # # from Askey "Ramanujan and Hypergeometric Series" in Berndt and # Rankin "Ramanujan: Essays and Surveys" # # this gives a sum of cosines of decreasing amp where the "k" # parameter determines the "index" (in FM nomenclature) -- higher k # = more cosines; the actual amount of the nth cos involves # hypergeometric series (looks like r^n/n! (~=e^n?) with a million # other terms). add_help(:kosine_summation, "kosine_summation(gen, r, k) \ Is a variant of sum-of-cosines; \ R controls successive sinusoid amplitude; \ K controls how many sinusoids are produced.") def kosine_summation(gen, r, k) r2 = r * r ((1.0 + r2) - 2 * r * oscil(gen)) ** -k * ((1.0 + r2) - 2 * r) ** k end alias make_kosine_summation make_oscil # legendre, fejer def fejer_sum(angle, n) if angle.zero? 1.0 else val = sin(0.5 * (n + 1) * angle) / (2.0 * sin(0.5 * angle)) 2.0 * ((val * val) / (n + 1)) end end def legendre_sum(angle, n) val = sin(angle * (n + 0.5)) / sin(0.5 * angle) val * val end # variations on sum-of-cosines # from "Trigonometric Delights" by Eli Maor def sum_of_n_sines(angle, n) a2 = angle * 0.5 den = sin(a2) if den.zero? 0.0 else (sin(n * a2) * sin((n + 1) * a2)) / den end end def sum_of_n_odd_sines(angle, n) angle = angle.to_f n = n.to_f den = sin(angle) na = sin(n * angle) if den.zero? 0.0 else (na * na) / den end end def sum_of_n_odd_cosines(angle, n) angle = angle.to_f n = n.to_f den = 2.0 * sin(angle) if den.zero? n else sin(2.0 * n * angle) / den end end # x = current phase, a = amp (more or less), N = 1..10 or # thereabouts, fi = phase increment Alexander Kritov suggests # time-varying "a" is good (this is a translation of his code) from # Stilson/Smith apparently -- was named "Discrete Summation Formula" # which doesn't convey anything to me def band_limited_sawtooth(x, a, n, fi) s4 = 1.0 + -2.0 * a * cos(x) + a * a if s4.zero? 0.0 else s1 = a ** (n - 1.0) * sin((n - 1.0) * x + fi) s2 = a ** n * sin(n * x + fi) s3 = a * sin(x + fi) (sin(fi) + -s3 + -s2 + s1) / s4 end end # brighten-slightly add_help(:brighten_slightly, "brighten_slightly(amount, snd=false, chn=false) \ Is a form of contrast-enhancement (AMOUNT between ca 0.1 and 1.0).") def brighten_slightly(amount, snd = false, chn = false) mx = maxamp brt = (TWO_PI * amount) / mx map_channel(lambda do |y| mx * sin(y * brt) end, 0, false, snd, chn, false, format("%s(%s", get_func_name, amount)) end add_help(:brighten_slightly_1, "brighten_slightly_1(coeffs, snd=false, chn=false) \ Is a form of contrast-enhancement: brighten_slightly-1([1, 0.5, 3, 1])") def brighten_slightly_1(coeffs, snd = false, chn = false) pcoeffs = partials2polynomial(coeffs) mx = maxamp(snd, chn) map_channel(lambda do |y| mx * polynomial(pcoeffs, y / mx) end, 0, false, snd, chn, false, format("%s(%s", get_func_name, coeffs)) end # FIR filters # Snd's (very simple) spectrum->coefficients procedure is: add_help(:spectrum2coeffs, "spectrum2coeffs(order, spectr) \ Returns FIR filter coefficients given the filter ORDER \ and desired spectral envelope (a vct).") def spectrum2coeffs(order, spectr) coeffs = make_vct(order) n = order m = ((n + 1) / 2.0).floor am = 0.5 * (n + 1) q = (PI * 2.0) / n jj = n - 1 m.times do |j| xt = 0.5 * spectr[0] (1...m).each do |i| xt = xt + spectr[i] * cos(q * i * (am - j - 1)) end coeff = 2.0 * (xt / n) coeffs[j] = coeff coeffs[jj] = coeff jj -= 1 end coeffs end add_help(:fltit_1, "fltit_1(order, spectrum) \ Creates an FIR filter from SPECTRUM and ORDER and \ returns a closure that calls it: \ map_channel(fltit_1(10, vct(0, 1.0, 0, 0, 0, 0, 0, 0, 1.0, 0)))") def fltit_1(order, spectr) flt = make_fir_filter(order, spectrum2coeffs(order, spectr)) lambda do |x| fir_filter(flt, x) end end # Hilbert transform add_help(:make_hilbert_transform, "make_hilbert_transform(len=30) \ Makes a Hilbert transform filter.") def make_hilbert_transform(len = 30) arrlen = len * 2 + 1 arr = make_vct(arrlen) (-len...len).each do |i| k = i + len denom = PI * i num = 1.0 - cos(denom) if num.zero? or i.zero? arr[k] = 0.0 else arr[k] = (num / denom) * (0.54 + 0.46 * cos(denom / len)) end end make_fir_filter(arrlen, arr) end add_help(:hilbert_transform, "hilbert_transform(f, in) \ Is the generator corresponding to make_hilbert_transform.") alias hilbert_transform fir_filter # highpass filter add_help(:make_highpass, "make_highpass(fc, len=30) \ Makes an FIR highpass filter.") def make_highpass(fc, len = 30) fc = fc.to_f arrlen = len * 2 + 1 arr = make_vct(arrlen) (-len...len).each do |i| k = i + len denom = PI * i num = -sin(fc * i) if i.zero? arr[k] = 1.0 - fc / PI else arr[k] = (num / denom) * (0.54 + 0.46 * cos((PI * i) / len)) end end make_fir_filter(arrlen, arr) end add_help(:highpass, "highpass(f, in) \ Is the generator corresponding to make_highpass.") alias highpass fir_filter # lowpass filter add_help(:make_lowpass, "make_lowpass(fc, len=30) \ Makes an FIR lowpass filter.") def make_lowpass(fc, len = 30) fc = fc.to_f arrlen = len * 2 + 1 arr = make_vct(arrlen) (-len...len).each do |i| k = i + len denom = PI * i num = sin(fc * i) if i.zero? arr[k] = fc / PI else arr[k] = (num / denom) * (0.54 + 0.46 * cos((PI * i) / len)) end end make_fir_filter(arrlen, arr) end add_help(:lowpass, "lowpass(f, in) \ Is the generator corresponding to make_lowpass.") alias lowpass fir_filter # bandpass filter add_help(:make_bandpass, "make_bandpass(flo, fhi, len=30) \ Makes an FIR bandpass filter.") def make_bandpass(flo, fhi, len = 30) flo = flo.to_f fhi = fhi.to_f arrlen = len * 2 + 1 arr = make_vct(arrlen) (-len...len).each do |i| k = i + len denom = PI * i num = sin(fhi * i) - sin(flo * i) if i.zero? arr[k] = (fhi - flo) / PI else arr[k] = (num / denom) * (0.54 + 0.46 * cos((PI * i) / len)) end end make_fir_filter(arrlen, arr) end add_help(:bandpass, "bandpass(f, in) \ Is the generator corresponding to make_bandpass.") alias bandpass fir_filter # bandstop filter add_help(:make_bandstop, "make_bandstop(flo, fhi, len=30) \ Makes an FIR bandstop (notch) filter.") def make_bandstop(flo, fhi, len = 30) flo = flo.to_f fhi = fhi.to_f arrlen = len * 2 + 1 arr = make_vct(arrlen) (-len...len).each do |i| k = i + len denom = PI * i num = sin(flo * i) - sin(fhi * i) if i.zero? arr[k] = 1.0 - (fhi - flo) / PI else arr[k] = (num / denom) * (0.54 + 0.46 * cos((PI * i) / len)) end end make_fir_filter(arrlen, arr) end add_help(:bandstop, "bandstop(f, in) \ Is the generator corresponding to make_bandstop.") alias bandstop fir_filter # differentiator add_help(:make_differentiator, "make_differentiator(len=30) \ Makes an FIR differentiator (highpass) filter.") def make_differentiator(len = 30) arrlen = len * 2 + 1 arr = make_vct(arrlen) (-len...len).each do |i| k = i + len if i.zero? arr[k] = 0.0 else arr[k] = (cos(PI * i) / i - sin(PI * i) / (PI * i * i)) * (0.54 + 0.46 * cos((PI * i) / len)) end end make_fir_filter(arrlen, arr) end add_help(:differentiator, "differentiator(f, in) \ Is the generator corresponding to make_differentiator.") alias differentiator fir_filter # IIR filters # # Butterworth filters (see also further below -- make-butter-lp et al) # # translated from CLM butterworth.cl: # # Sam Heisz, January 1998 # inspired by some unit generators written for Csound by Paris # Smaragdis who based his work on formulas from Charles Dodge, # Computer music: synthesis, composition, and performance. add_help(:butter, "butter(b, sig) \ Is the generator side for the various make_butter procedure.") alias butter filter add_help(:make_butter_high_pass, "make_butter_high_pass(freq) \ Makes a Butterworth filter with high pass cutoff at FREQ.") def make_butter_high_pass(freq) r = tan(PI * freq / srate()) r2 = r * r c1 = 1.0 / (1.0 + r * sqrt(2.0) + r2) c2 = -2.0 * c1 c3 = c1 c4 = 2.0 * (r2 - 1.0) * c1 c5 = ((1.0 - r * sqrt(2.0)) + r2) * c1 make_filter(3, vct(c1, c2, c3), vct(0.0, c4, c5)) end add_help(:make_butter_low_pass, "make_butter_low_pass(freq) \ Makes a Butterworth filter with high pass cutoff at FREQ. \ The result can be used directly: \ filter_sound(make_butter_low_pass(500.0)), or via the 'butter' generator.") def make_butter_low_pass(freq) r = 1.0 / tan(PI * freq / srate()) r2 = r * r c1 = 1.0 / (1.0 + r * sqrt(2.0) + r2) c2 = 2.0 * c1 c3 = c1 c4 = 2.0 * (1.0 - r2) * c1 c5 = ((1.0 - r * sqrt(2.0)) + r2) * c1 make_filter(3, vct(c1, c2, c3), vct(0.0, c4, c5)) end add_help(:make_butter_band_pass, "make_butter_band_pass(freq, band) \ Makes a bandpass Butterworth filter with low edge at FREQ and width BAND.") def make_butter_band_pass(freq, band) d = 2.0 * cos(2.0 * PI * freq / srate()) c = 1.0 / tan(PI * band / srate()) c1 = 1.0 / (1.0 + c) c2 = 0.0 c3 = -c1 c4 = -c * d * c1 c5 = (c - 1.0) * c1 make_filter(3, vct(c1, c2, c3), vct(0.0, c4, c5)) end add_help(:make_butter_band_reject, "make_butter_band_reject(freq, band) \ Makes a band-reject Butterworth filter with low edge at FREQ and width BAND.") def make_butter_band_reject(freq, band) d = 2.0 * cos(2.0 * PI * freq / srate()) c = tan(PI * band / srate()) c1 = 1.0 / (1.0 + c) c2 = -d * c1 c3 = c1 c4 = c2 c5 = (1.0 - c) * c1 make_filter(3, vct(c1, c2, c3), vct(0.0, c4, c5)) end # from "DSP Filter Cookbook" by Lane et al, Prompt Pubs, 2001 # # use with the filter generator # (define gen (make-iir-high-pass-2 1000)) # (filter gen 1.0) # etc add_help(:make_biquad, "make_biquad(a0, a1, a2, b1, b2) \ Returns a biquad filter (use with the CLM filter gen).") def make_biquad(a0, a1, a2, b1, b2) make_filter(3, vct(a0, a1, a2), vct(0.0, b1, b2)) end # din=(sqrt 2.0) for example (suggested range 0.2...10) def make_iir_low_pass_2(fc, din = false) fc = fc.to_f theta = (TWO_PI * fc) / mus_srate() d = (din or sqrt(2.0)) beta = 0.5 * ((1.0 - (d / 2.0) * sin(theta)) / (1.0 + (d / 2.0) * sin(theta))) gamma = (0.5 + beta) * cos(theta) alpha = 0.5 * (0.5 + beta + -gamma) make_filter(3, vct(alpha, 2.0 * alpha, alpha), vct(0.0, -2.0 * gamma, 2.0 * beta)) end def make_iir_high_pass_2(fc, din = false) fc = fc.to_f theta = (TWO_PI * fc) / mus_srate() d = (din or sqrt(2.0)) beta = 0.5 * ((1.0 - (d / 2.0) * sin(theta)) / (1.0 + (d / 2.0) * sin(theta))) gamma = (0.5 + beta) * cos(theta) alpha = 0.5 * (0.5 + beta + gamma) make_filter(3, vct(alpha, -2.0 * alpha, alpha), vct(0.0, -2.0 * gamma, 2.0 * beta)) end def make_iir_band_pass_2(f1, f2) f1 = f1.to_f f2 = f2.to_f theta = (TWO_PI * sqrt(f1 * f2)) / mus_srate() q = sqrt(f1 * f2) / (f2 - f1) t2 = tan(theta / (2 * q)) beta = 0.5 * ((1.0 - t2) / (1.0 + t2)) gamma = (0.5 + beta) * cos(theta) alpha = 0.5 - beta make_filter(3, vct(alpha, 0.0, -alpha), vct(0.0, -2.0 * gamma, 2.0 * beta)) end def make_iir_band_stop_2(f1, f2) f1 = f1.to_f f2 = f2.to_f theta = (TWO_PI * sqrt(f1 * f2)) / mus_srate() q = sqrt(f1 * f2) / (f2 - f1) t2 = tan(theta / (2 * q)) beta = 0.5 * ((1.0 - t2) / (1.0 + t2)) gamma = (0.5 + beta) * cos(theta) alpha = 0.5 + beta make_filter(3, vct(alpha, -2.0 * gamma, alpha), vct(0.0, -2.0 * gamma, 2.0 * beta)) end def make_eliminate_hum(hum_freq = 60.0, hum_harmonics = 5, bandwidth = 10) b2 = 0.5 * bandwidth make_array(hum_harmonics) do |i| center = (i + 1.0) * hum_freq make_iir_band_stop_2(center - b2, center + b2) end end def eliminate_hum(gens, x0) val = x0 gens.each do |gen| val = filter(gen, val) end val end # bandpass, m is gain at center of peak # use map-channel with this one (not clm-channel or filter) def make_peaking_2(f1, f2, m) f1 = f1.to_f f2 = f2.to_f theta = (TWO_PI * sqrt(f1 * f2)) / mus_srate() q = sqrt(f1 * f2) / (f2 - f1) t2 = (4.0 / (m + 1.0)) * tan(theta / (2 * q)) beta = 0.5 * ((1.0 - t2) / (1.0 + t2)) gamma = (0.5 + beta) * cos(theta) alpha = 0.5 - beta flt = make_filter(3, vct(alpha, 0.0, -alpha), vct(0.0, -2.0 * gamma, 2.0 * beta)) lambda do |x| x + (m - 1.0) * filter(flt, x) end end # convert cascade coeffs to canonical form # from Orfanidis "Introduction to Signal Processing def c2c_conv(m, h, l, x, y) (l + m).times do |i| y[i] = 0.0 ([0, i - (1 + l)].max..[i, m].min).each do |j| y[i] = y[i] + h[j] * x[i - j] end end end add_help(:cascade2canonical, "cascade2canonical(a) \ Converts a list of cascade coeffs (vcts with 3 entries) to canonical form.") def cascade2canonical(a) k = a.length d = make_vct(2 * k + 1) a1 = make_vct(2 * k + 1) a1[0] = 1.0 k.times do |i| c2c_conv(2, a[i], 2 * i + 1, a1, d) (2 * i + 3).times do |j| a1[j] = d[j] end end a1 end # order is M*2, fc is cutoff freq (Hz) add_help(:make_butter_lp, "make_butter_lp(m, fc) \ Returns a butterworth low-pass filter; \ its order is M * 2, FC is the cutoff frequency in Hz.") def make_butter_lp(m, fc) fc = fc.to_f xcoeffs = make_array(m) ycoeffs = make_array(m) theta = (TWO_PI * fc) / mus_srate() st = sin(theta) ct = cos(theta) m.times do |k| d = 2.0 * sin((PI * (2.0 * (k + 1.0) - 1.0)) / (4.0 * m)) beta = 0.5 * ((1.0 - 0.5 * d * st) / (1.0 + 0.5 * d * st)) gamma = ct * (0.5 + beta) alpha = 0.25 * (0.5 + beta + -gamma) xcoeffs[k] = vct(2.0 * alpha, 4.0 * alpha, 2.0 * alpha) ycoeffs[k] = vct(1.0, -2.0 * gamma, 2.0 * beta) end make_filter(2 * m + 1, cascade2canonical(xcoeffs), cascade2canonical(ycoeffs)) end # order is M*2, fc is cutoff freq (Hz) add_help(:make_butter_hp, "make_butter_hp(m, fc) \ Returns a butterworth high-pass filter; \ its order is M * 2, FC is the cutoff frequency in Hz.") def make_butter_hp(m, fc) fc = fc.to_f xcoeffs = make_array(m) ycoeffs = make_array(m) theta = (TWO_PI * fc) / mus_srate() st = sin(theta) ct = cos(theta) m.times do |k| d = 2.0 * sin((PI * (2.0 * (k + 1.0) - 1.0)) / (4.0 * m)) beta = 0.5 * ((1.0 - 0.5 * d * st) / (1.0 + 0.5 * d * st)) gamma = ct * (0.5 + beta) alpha = 0.25 * (0.5 + beta + gamma) xcoeffs[k] = vct(2.0 * alpha, -4 * alpha, 2.0 * alpha) ycoeffs[k] = vct(1.0, -2.0 * gamma, 2.0 * beta) end make_filter(2 * m + 1, cascade2canonical(xcoeffs), cascade2canonical(ycoeffs)) end # order is M*2, f1 and f2 are band edge freqs (Hz) add_help(:make_butter_bp, "make_butter_bp(m, f1, f2) \ Returns a butterworth band-pass filter; \ its order is M * 2, F1 and F2 are the band edge frequencies in Hz.") def make_butter_bp(m, f1, f2) f1 = f1.to_f f2 = f2.to_f xcoeffs = make_array(m) ycoeffs = make_array(m) f0 = sqrt(f1 * f2) q = f0 / (f2 - f1) theta0 = (TWO_PI * f0) / mus_srate() de = (2.0 * tan(theta0 / (2.0 * q))) / sin(theta0) de2 = de / 2.0 tn0 = tan(theta0 * 0.5) k = j = 1 m.times do |i| dk = 2.0 * sin((PI * (2.0 * k - 1.0)) / (2.0 * m)) ak = (1.0 + de2 * de2) / (dk * de2) dk1 = sqrt((de * dk) / (ak + sqrt(ak * ak - 1.0))) bk = de2 * (dk / dk1) wk = (bk + sqrt(bk * bk - 1.0)).real thetajk = ((j == 1) ? (2.0 * atan(tn0 / wk)) : (2.0 * atan(tn0 * wk))) betajk = 0.5 * ((1.0 - 0.5 * dk1 * sin(thetajk)) / (1.0 + 0.5 * dk1 * sin(thetajk))) gammajk = (0.5 + betajk) * cos(thetajk) wk2 = (wk - 1.0 / wk) / dk1 alphajk = 0.5 * (0.5 - betajk) * sqrt(1.0 + wk2 * wk2) xcoeffs[i] = vct(2.0 * alphajk, 0.0, -2.0 * alphajk) ycoeffs[i] = vct(1.0, -2.0 * gammajk, 2.0 * betajk) if j == 1 j = 2 else k += 1 j = 1 end end make_filter(2 * m + 1, cascade2canonical(xcoeffs), cascade2canonical(ycoeffs)) end # order is M*2, f1 and f2 are band edge freqs (Hz) add_help(:make_butter_bs, "make_butter_bs(m, f1, f2) \ Returns a butterworth band-stop filter; \ its order is M * 2, F1 and F2 are the band edge frequencies in Hz.") def make_butter_bs(m, f1, f2) f1 = f1.to_f f2 = f2.to_f xcoeffs = make_array(m) ycoeffs = make_array(m) f0 = sqrt(f1 * f2) q = f0 / (f2 - f1) theta0 = (TWO_PI * f0) / mus_srate() de = (2.0 * tan(theta0 / (2.0 * q))) / sin(theta0) de2 = de / 2.0 ct = cos(theta0) tn0 = tan(theta0 * 0.5) k = j = 1 m.times do |i| dk = 2.0 * sin((PI * (2.0 * k - 1.0)) / (2.0 * m)) ak = (1.0 + de2 * de2) / (dk * de2) dk1 = sqrt((de * dk) / (ak + sqrt(ak * ak - 1.0))) bk = de2 * (dk / dk1) wk = (bk + sqrt(bk * bk - 1.0)).real thetajk = ((j == 1) ? (2.0 * atan(tn0 / wk)) : (2.0 * atan(tn0 * wk))) betajk = 0.5 * ((1.0 - 0.5 * dk1 * sin(thetajk)) / (1.0 + 0.5 * dk1 * sin(thetajk))) gammajk = (0.5 + betajk) * cos(thetajk) alphajk = 0.5 * (0.5 + betajk) * ((1.0 - cos(thetajk)) / (1.0 - ct)) xcoeffs[i] = vct(2.0 * alphajk, -4.0 * ct * alphajk, 2.0 * alphajk) ycoeffs[i] = vct(1.0, -2.0 * gammajk, 2.0 * betajk) if j == 1 j = 2 else k += 1 j = 1 end end make_filter(2 * m + 1, cascade2canonical(xcoeffs), cascade2canonical(ycoeffs)) end # notch filters def make_notch_frequency_response(cur_srate, freqs, notch_width = 2) cur_srate = cur_srate.to_f notch_width = notch_width.to_f freq_response = [0.0, 1.0] freqs.each do |f| # left upper y hz freq_response.push((2.0 * (f - notch_width)) / cur_srate) # left upper y resp freq_response.push(1.0) # left bottom y hz freq_response.push((2.0 * (f - notch_width / 2.0)) / cur_srate) # left bottom y resp freq_response.push(0.0) # right bottom y hz freq_response.push((2.0 * (f + notch_width / 2.0)) / cur_srate) # right bottom y resp freq_response.push(0.0) # right upper y hz freq_response.push((2.0 * (f + notch_width)) / cur_srate) # right upper y resp freq_response.push(1.0) end freq_response.push(1.0, 1.0) end add_help(:notch_channel, "notch_channel(freqs, order=false, beg=false, dur=false, \ snd=false, chn=false, edpos=false, trunc=true, notch_width=2) \ Returns notch filter removing freqs.") def notch_channel(freqs, filter_order = false, beg = false, dur = false, snd = false, chn = false, edpos = false, truncate = true, notch_width = 2) sr = srate(snd).to_f lm = [framples(snd, chn), 2 ** (log(sr / notch_width) / log(2.0)).floor].min filter_channel(make_notch_frequency_response(sr, freqs, notch_width), (filter_order or lm), beg, dur, snd, chn, edpos, truncate, format("%s(%p, %s, %s, %s", get_func_name, freqs, filter_order, beg, dur)) end add_help(:notch_sound, "notch_sound(freqs, order=false, \ snd=false, chn=false, notch_width=2) \ Returns notch filter removing freqs.") def notch_sound(freqs, filter_order = false, snd = false, chn = false, notch_width = 2) sr = srate(snd).to_f lm = [framples(snd, chn), 2 ** (log(sr / notch_width) / log(2.0)).floor].min filter_sound(make_notch_frequency_response(sr, freqs, notch_width), (filter_order or lm), snd, chn, false, format("%s(%p, %s, 0, false", get_func_name, freqs, filter_order)) end add_help(:notch_selection, "notch_selection(freqs, order=false, notch_width=2) \ Returns notch filter removing freqs.") def notch_selection(freqs, filter_order = false, snd = false, chn = false, notch_width = 2) if selection? sr = selection_srate.to_f fr = selection_framples() lm = [fr, 2 ** (log(sr / notch_width) / log(2.0)).floor].min filter_selection(make_notch_frequency_response(sr, freqs, notch_width), (filter_order or lm)) end end # fractional Fourier Transform, z transform # # translated from the fxt package of Joerg Arndt add_help(:fractional_fourier_transform, "fractional_fourier_transform(real, imaginary, n, angle) \ Performs a fractional Fourier transform on data; \ if angle=1.0, you get a normal Fourier transform.") def fractional_fourier_transform(fr, fi, n, v) hr = make_vct(n) hi = make_vct(n) ph0 = (v * TWO_PI) / n n.times do |w| sr = 0.0 si = 0.0 n.times do |k| phase = ph0 * k * w c = cos(phase) s = sin(phase) x = fr[k] y = fi[k] r = x * c - y * s i = y * c + x * s sr += r si += i hr[w] = sr hi[w] = si end end [hr, hi] end # using vector to allow complex sums (z=e^2*pi*i/n -> fourier transform) # z_transform(data, n, exp(Complex(0.0, (2.0 / n) * PI))) add_help(:z_transform, "z_transform(data, n, z) \ Performs a Z transform on data; \ if z=e^2*pi*j/n you get a Fourier transform; \ complex results in returned vector.") def z_transform(f, n, z) make_array(n) do |w| sum = 0.0 t = 1.0 m = z ** w n.times do |k| sum = sum + f[k] * t t *= m end sum end end # slow Hartley transform # # taken from Perry Cook's SignalProcessor.m (the slow version of the # Hartley transform) add_help(:dht, "dht(data) \ Returns the Hartley transform of DATA.") def dht(data) len = data.length arr = make_vct(len) w = TWO_PI / len len.times do |i| data.each_with_index do |val, j| arr[i] = arr[i] + val * (cos(i * j * w) + sin(i * j * w)) end end arr end add_help(:find_sine, "find_sine(freq, beg, dur, snd=false) \ Returns the amplitude and initial-phase (for sin) at FREQ between BEG and DUR.") def find_sine(freq, beg, dur, snd = false) incr = (TWO_PI * freq) / srate(snd) sw = 0.0 cw = 0.0 reader = make_sampler(beg, snd) dur.times do |i| samp = next_sample(reader) inc = i * incr sw = sw + samp * sin(inc) cw = cw + samp * cos(inc) end [2.0 * (sqrt(sw * sw + cw * cw) / dur), atan2(cw, sw)] end # this is a faster version of find-sine using the "Goertzel algorithm" # taken from R Lyons "Understanding DSP" p 529 # it returns the same result as find_sine above if you take (* 2 (/ # (goertzel...) dur)) -- see snd-test.rb examples add_help(:goertzel, "goertzel(freq, beg=0, dur=false, snd=false) \ Returns the amplitude of the FREQ spectral component.") def goertzel(freq, beg = 0, dur = false, snd = false) y0 = 0.0 y1 = 0.0 y2 = 0.0 rfreq = (TWO_PI * freq) / srate(snd) cs = 2.0 * cos(rfreq) scan_channel(lambda do |y| y2, y1 = y1, y0 y0 = (y1 * cs - y2) + y false end, beg, (dur or framples(snd)), snd) (y0 - y1 * exp(Complex(0.0, -rfreq))).abs end add_help(:make_spencer_filter, "make_spencer_filter() \ Is a version of make_fir_filter; \ it returns one of the standard smoothing filters from \ the era when computers were human beings.") def make_spencer_filter data = vct(-3, -6, -5, 3, 21, 46, 67, 74, 67, 46, 21, 3, -5, -6, -3) data.map! do |n| n / 320.0 end make_fir_filter(data.length, data) end # any-random # # arbitrary random number distributions via the "rejection method" def any_random(amount, e = false) if amount.zero? 0.0 else unless e random(amount) else next_random = lambda do | | len = e.length x = random(e[len - 2].to_f) y = random(1.0) if y <= envelope_interp(x, e) x else next_random.call end end.call end end end def gaussian_distribution(s) e = [] den = 2.0 * s * s x = 0.0 y = -4.0 21.times do |i| e.push(exp(-((y * y) / den))) e.push(x) x += 0.05 y += 0.4 end e end def pareto_distribution(a) e = [] scl = 1.0 ** (a + 1.0) / a x = 0.0 y = 1.0 21.times do |i| e.push(scl * (a / y ** (a + 1.0))) e.push(x) x += 0.05 y += 0.2 end e end # uniform distribution # map_channel(lambda do |y| any_random(1.0, [0, 1, 1, 1])) # mostly toward 1.0 # map_channel(lambda do |y| any_random(1.0, [0, 0, 0.95, 0.1, 1, 1])) # let(gaussian-distribution(1.0)) do |g| # map_channel(lambda do |y| any_random(1.0, g)) # end # let(pareto-distribution(1.0)) do |g| # map_channel(lambda do |y| any_random(1.0, g)) # end # this is the inverse integration function used by CLM to turn a # distribution function into a weighting function def inverse_integrate(dist, data_size = 512, e_size = 50) first_sum = sum = dist[1].to_f x0 = dist[0].to_f x1 = dist[-2].to_f xincr = (x1 - x0) / e_size.to_f x = x0 e = make_array(e_size * 2) 0.step(e_size * 2 - 1, 2) do |i| e[i + 1] = x e[i] = sum sum += envelope_interp(x, dist) x += xincr end incr = (e[-2] - first_sum) / (data_size - 1) x = first_sum - incr make_vct!(data_size) do x += incr envelope_interp(x, e) end end def gaussian_envelope(s) den = 2.0 * s * s x = -1.0 y = -4.0 e = make_array(42) 0.step(41, 2) do |i| e[i] = x e[i + 1] = exp(-((y * y) / den)) x += 0.1 y += 0.4 end e end # make_rand(:envelope, gaussian-envelope(1.0)) # Julius Smith stuff # # these are from "Mathematics of the DFT", W3K Pubs add_help(:channel_mean, "channel_mean(snd, chn) \ Returns the average of the samples in the given channel: /n") def channel_mean(snd, chn) sum = 0.0 n = framples(snd, chn) scan_channel(lambda do |y| sum += y false end, 0, n, snd, chn) sum / n end add_help(:channel_total_energy, "channel_total_energy(snd, chn) \ Returns the sum of the squares of all the samples in the given channel: ") def channel_total_energy(snd, chn) sum = 0.0 scan_channel(lambda do |y| sum = sum + y * y false end, 0, framples(snd, chn), snd, chn) sum end add_help(:channel_average_power, "channel_average_power(snd, chn) \ Returns the average power in the given channel: /n") def channel_average_power(snd, chn) channel_total_energy(snd, chn) / framples(snd, chn) end add_help(:channel_rms, "channel_rms(snd, chn) \ Returns the RMS value of the samples in the given channel: sqrt(/n)") def channel_rms(snd, chn) sqrt(channel_average_power(snd, chn)) end add_help(:channel_variance, "channel_variance(snd, chn) \ Returns the sample variance in the given channel: -((/ n)^2") def channel_variance(snd, chn) n = framples(snd, chn).to_f mu = (n / (n - 1.0)) * channel_mean(snd, chn) p = channel_total_energy(snd, chn) p - mu * mu end add_help(:channel_norm, "channel_norm(snd, chn) \ Returns the norm of the samples in the given channel: sqrt()") def channel_norm(snd, chn) sqrt(channel_total_energy(snd, chn)) end add_help(:channel_lp, "channel_lp(p, snd, chn) \ Returns the Lp norm of the samples in the given channel.") def channel_lp(lp, snd, chn) sum = 0.0 scan_channel(lambda do |y| sum = sum + y.abs ** lp false end, 0, framples(snd, chn), snd, chn) sum ** (1.0 / lp) end add_help(:channel_lp_inf, "channel_lp_inf(snd, chn) \ Returns the maxamp in the given channel (the name is \ just math jargon for maxamp).") def channel_lp_inf(snd, chn) mx = 0.0 scan_channel(lambda do |y| mx = [mx, y.abs].max false end, 0, framples(snd, chn), snd, chn) mx end add_help(:channel2_inner_product, "channel2_inner_product(s1, c1, s2, c2) \ Returns the inner-product of the two channels: ") def channel2_inner_product(s1, c1, s2, c2) sum = 0.0 r1 = make_sampler(0, s1, c1) r2 = make_sampler(0, s2, c2) framples(s1, c1).times do |i| sum = sum + r1.call * r2.call end sum end add_help(:channel2_angle, "channel2_angle(s1, c1, s2, c2) \ Treats the two channels as vectors, \ returning the ANGLE between them: acos(/(sqrt()*sqrt()))") def channel2_angle(s1, c1, s2, c2) inprod = channel2_inner_product(s1, c1, s2, c2) norm1 = channel_norm(s1, c1) norm2 = channel_norm(s2, c2) acos(inprod / (norm1 * norm2)) end if defined? acos add_help(:channel2_orthogonal?, "channel2_orthogonal?(s1, c1, s2, c2) \ Returns true if the two channels' inner-product is 0: ==0") def channel2_orthogonal?(s1, c1, s2, c2) channel2_inner_product(s1, c1, s2, c2).zero? end add_help(:channel2_coefficient_of_projection, "channel2_coefficient_of_projection(s1, c1, s2, c2) \ Returns /") def channel2_coefficient_of_projection(s1, c1, s2, c2) channel2_inner_product(s1, c1, s2, c2) / channel_total_energy(s1, c1) end # end of JOS stuff add_help(:channel_distance, "channel_distance(s1, c1, s2, c2) \ Returns the euclidean distance between the two channels: sqrt()") def channel_distance(s1, c1, s2, c2) r1 = make_sampler(0, s1, c1) r2 = make_sampler(0, s2, c2) sum = 0.0 [framples(s1, c1), framples(s2, c2)].min.times do diff = r1.call - r2.call sum = sum + diff * diff end sqrt(sum) end add_help(:periodogram, "periodogram(n) \ Displays an N point Bartlett periodogram of the samples in the current channel") def periodogram(n) len = framples() average_data = make_vct(n) rd = make_sampler(0) n2 = n * 2 rl = make_vct(n2) im = make_vct(n2) len.times do vct_scale!(rl, 0.0) vct_scale!(im, 0.0) n.times do |k| rl[k] = rd.call end mus_fft(rl, im) n.times do |k| average_data[k] = average_data[k] + rl[k] * rl[k] + im[k] * im[k] end end graph(vct_scale!(average_data, 1.0 / (len.to_f / n).ceil)) end # ssb-am friends add_help(:shift_channel_pitch, "shift_channel_pitch(freq, order=40, beg=0, dur=false, \ snd=false, chn=false, edpos=false) \ Uses the ssb-am CLM generator to shift the given channel \ in pitch without changing its length. \ The higher ORDER, the better usually.") def shift_channel_pitch(freq, order = 40, beg = 0, dur = false, snd = false, chn = false, edpos = false) gen = make_ssb_am(freq, order) map_channel(lambda do |y| ssb_am(gen, y) end, beg, dur, snd, chn, edpos, format("%s(%s, %s, %s, %s", get_func_name, freq, order, beg, dur)) end add_help(:hz_to_2pi, "hz_to_2pi(freq) \ Is like hz2radians but uses the current sound's srate, not mus_srate.") def hz_to_2pi(freq) (TWO_PI * freq) / srate() end def ssb_bank(old_freq, new_freq, pairs, order = 40, bw = 50.0, beg = 0, dur = false, snd = false, chn = false, edpos = false) factor = (new_freq - old_freq.to_f) / old_freq mx = maxamp ssbs = make_array(pairs) bands = make_array(pairs) do |i| aff = (i + 1.0) * old_freq bwf = bw * (1.0 + (i + 1.0) / (2.0 * pairs)) ssbs[i] = make_ssb_am((i + 1.0) * factor * old_freq) make_bandpass(hz_to_2pi(aff - bwf), hz_to_2pi(aff + bwf), order) end as_one_edit_rb("%s(%s, %s, %s, %s, %s, %s, %s", get_func_name, old_freq, new_freq, pairs, order, bw, beg, dur) do | | nmx = 0.0 map_channel_rb(beg, dur, snd, chn, edpos) do |y| sum = 0.0 ssbs.zip(bands) do |sbs, bds| sum += ssb_am(sbs, bandpass(bds, y)) end nmx = [nmx, sum.abs].max sum end scale_channel(mx / nmx, beg, dur, snd, chn) end end # this version adds a frequency envelope # ssb_bank_env(557, 880, [0, 0, 1, 100.0], 7) def ssb_bank_env(old_freq, new_freq, freq_env, pairs, order = 40, bw = 50.0, beg = 0, dur = false, snd = false, chn = false, edpos = false) factor = (new_freq - old_freq.to_f) / old_freq mx = maxamp ssbs = make_array(pairs) frenvs = make_array(pairs) bands = make_array(pairs) do |i| aff = (i + 1.0) * old_freq bwf = bw * (1.0 + (i + 1.0) / (2.0 * pairs)) ssbs[i] = make_ssb_am((i + 1.0) * factor * old_freq) frenvs[i] = make_env(:envelope, freq_env, :scaler, hz2radians(i.to_f), :length, framples(snd, chn) - 1) make_bandpass(hz_to_2pi(aff - bwf), hz_to_2pi(aff + bwf), order) end as_one_edit_rb("%s(%s, %s, %s, %s, %s, %s, %s, %s", get_func_name, old_freq, new_freq, freq_env.inspect, pairs, order, bw, beg, dur) do | | nmx = 0.0 map_channel_rb(beg, dur, snd, chn, edpos) do |y| sum = 0.0 ssbs.each_with_index do |sbs, i| sum += ssb_am(sbs, bandpass(bands[i], y), env(frenvs[i])) end nmx = [nmx, sum.abs].max sum end scale_channel(mx / nmx, beg, dur, snd, chn) end end #vct|channel|spectral-polynomial def vct_polynomial(v, coeffs) new_v = Vct.new(v.length, coeffs.last) (coeffs.length - 2).downto(0) do |i| new_v.multiply!(v).offset!(coeffs[i]) end new_v end def channel_polynomial(coeffs, snd = false, chn = false) len = framples(snd, chn) v = channel2vct(0, len, snd, chn) vct2channel(vct_polynomial(v, coeffs), 0, len, snd, chn, false, format("%s(%s", get_func_name, coeffs.to_str)) end # channel_polynomial(vct(0.0, 0.5)) = x*0.5 # channel_polynomial(vct(0.0, 1.0, 1.0, 1.0)) = x*x*x + x*x + x # # convolution -> * in freq def spectral_polynomial(coeffs, snd = false, chn = false) len = framples(snd, chn) sound = channel2vct(0, len, snd, chn) num_coeffs = coeffs.length fft_len = if num_coeffs < 2 len else (2.0 ** (log((num_coeffs - 1.0) * len) / log(2.0)).ceil).to_i end rl1 = make_vct(fft_len) rl2 = make_vct(fft_len) new_sound = make_vct(fft_len) if coeffs[0] > 0.0 new_sound.map! do mus_random(coeffs[0]) end end if num_coeffs > 1 new_sound.add!(sound.scale(coeffs[1])) if num_coeffs > 2 peak = maxamp(snd, chn) rl1.scale!(0.0).add!(sound) (2...num_coeffs).each do |i| convolution(rl1, rl2.scale!(0.0).add(sound), fft_len) new_sound.add!(rl1.scale((coeffs[i] * peak) / rl1.peak)) end new_sound.scale!(peak / new_sound.peak) end end vct2channel(new_sound, 0, [len, len * (num_coeffs - 1)].max, snd, chn, false, format("%s(%s", get_func_name, coeffs.to_str)) end # SCENTROID # # by Bret Battey # Version 1.0 July 13, 2002 # translated to Snd/Scheme Bill S 19-Jan-05 # # Returns the continuous spectral centroid envelope of a sound. The # spectral centroid is the "center of gravity" of the spectrum, and it # has a rough correlation to our sense of "brightness" of a sound. # # [Beauchamp, J., "Synthesis by spectral amplitude and 'brightness' matching # analyzed musical sounds". # Journal of Audio Engineering Society 30(6), 396-406] # # The formula used is: # C = [SUMF(n)A(n)] / [SUMA(n)] # Where j is the number of bins in the analysis, # F(n) is the frequency of a given bin, # A(n) is the magnitude of the given bin. # # If a pitch envelope for the analyzed sound is available, the results # of SCENTROID can be used with the function NORMALIZE-CENTROID, # below, to provide a "normalized spectral centroid". # # DB-FLOOR -- Frames below this decibel level (0 dB = max) will be # discarded and returned with spectral centroid = 0 # # RFREQ -- Rendering frequency. Number of measurements per second. # # FFTSIZE -- FFT window size. Must be a power of 2. 4096 is # recommended. add_help(:scentroid, "scentroid(file, beg=0, dur=false, db_floor=-40, \ rfreq=100, fftsize=4096) \ Returns the spectral centroid envelope of a sound; \ RFREQ is the rendering frequency, the number of measurements per second; \ DB_FLOOR is the level below which data will be ignored.") def scentroid(file, beg = 0.0, dur = false, db_floor = -40, rfreq = 100.0, fftsize = 4096) assert_type(File.exist?(file), file, 0, "an existing file") fsr = srate(file) incrsamps = (fsr / rfreq).floor start = (beg * fsr).floor ende = (start + (dur ? (dur * fsr).floor : (framples(file) - beg))).floor fdr = make_vct(fftsize) fdi = make_vct(fftsize) windows = ((ende - start) / incrsamps).floor + 1 results = make_vct(windows) fft2 = (fftsize / 2.0).floor binwidth = fsr / fftsize.to_f rd = make_readin(file) loc = 0 start.step(ende, incrsamps) do |i| rd.location = i sum_of_squares = 0.0 fdr.map! do val = readin(rd) sum_of_squares = sum_of_squares + val * val val end if linear2db(sqrt(sum_of_squares / fftsize.to_f)) >= db_floor numsum = 0.0 densum = 0.0 fdi.map! do |x| 0.0 end mus_fft(fdr, fdi, fftsize) rectangular2polar(fdr, fdi) fft2.times do |j| numsum = numsum + j * binwidth * fdr[j] densum += fdr[j] end results[loc] = numsum / densum end loc += 1 end results end # invert_filter inverts an FIR filter # # say we previously filtered a sound via filter_channel(vct(0.5, 0.25, 0.125)) # and we want to undo it without using undo_edit: # filter_channel(invert_filter(vct(0.5, 0.25, 0.125))) # # there are a million gotchas here. The primary one is that the inverse # filtercan "explode" -- the coefficients can grow without bound. For # example, any filter returned by spectrum2coeffs above will be a problem # (it always returns a "linear phase" filter). add_help(:invert_filter, "invert_filter(coeffs) \ Tries to return an inverse filter to undo the effect of the FIR filter coeffs.") def invert_filter(fcoeffs) order = fcoeffs.length + 32 coeffs = Vct.new(order) fcoeffs.each_with_index do |val, i| coeffs[i] = val end nfilt = Vct.new(order) nfilt[0] = 1.0 / coeffs.first (1...order).each do |i| sum = 0.0 k = i i.times do |j| sum = sum + nfilt[j] * coeffs[k] k -= 1 nfilt[i] = sum / -coeffs.first end end nfilt end # Volterra filter # # one of the standard non-linear filters # this version is taken from Monson Hayes "Statistical DSP and Modeling" # it is a slight specialization of the form mentioned by J O Smith and others class Volterra_filter < Musgen def initialize(acoeffs, bcoeffs) super() @as = acoeffs @bs = bcoeffs @xs = Vct.new([acoeffs.length, bcoeffs.length].max) end def inspect format("%s.new(%s, %s)", self.class, @as.to_str, @bs.to_str) end def to_s format("#<%s acoeffs: %s, bcoeffs: %s>", self.class, @as, @bs) end def run_func(val1 = 0.0, val2 = 0.0) volterra_filter(val1) end def volterra_filter(x) xlen = @xs.length @xs.move!(xlen - 1, xlen - 2, true) @xs.first = x sum = dot_product(@as, @xs, @as.length) @bs.length.times do |i| @bs.length.times do |j| sum = sum + @bs[j] * @xs[i] * @xs[j] end end sum end end add_help(:make_volterra_filter, "make_volterra_filter(acoeffs, bcoeffs) \ Returns a list for use with volterra-filter, \ producing one of the standard non-linear filters.") def make_volterra_filter(acoeffs, bcoeffs) Volterra_filter.new(acoeffs, bcoeffs) end add_help(:volterra_filter, "volterra_filter(flt, x) \ Takes FLT, a Volterra_filter object returned by make_volterra_filter, \ and an input X, and returns the (non-linear filtered) result.") def volterra_filter(flt, x) flt.volterra_filter(x) end # flt = make_volterra_filter(vct(0.5, 0.1), vct(0.3, 0.2, 0.1)) # map_channel(lambda do |y| volterra_filter(flt, y) end) # moving-sum generator (the sum norm or 1-norm) add_help(:make_moving_sum, "make_moving_sum(size=128) \ Returns a moving-sum generator. \ The generator keeps a running window of the last SIZE inputs, \ returning the sum of the absolute values of the samples in that window.") def make_moving_sum(size = 128) gen = make_moving_average(size) gen.increment = 1.0 gen end add_help(:moving_sum, "moving_sum(gen, y) \ Returns the sum of the absolute values in a moving \ window over the last few inputs.") def moving_sum(gen, y) moving_average(gen, y.abs) end def make_unmoving_sum() make_one_pole(1.0, -1.0) end alias unmoving_sum one_pole # moving-rms generator add_help(:make_moving_rms, "make_moving_rms(size=128) \ Returns a moving-rms generator. \ The generator keeps a running window of the last SIZE inputs, \ returning the rms of the samples in that window.") def make_moving_rms(size = 128) make_moving_average(size) end add_help(:moving_rms, "moving_rms(gen, y) \ Returns the rms of the values in a window over the last few inputs.") def moving_rms(gen, y) sqrt([0.0, moving_average(gen, y * y)].max) end # moving-length generator (euclidean norm or 2-norm) add_help(:make_moving_length, "make_moving_length(size=128) \ Returns a moving-length generator. \ The generator keeps a running window of the last SIZE inputs, \ returning the euclidean length of the vector in that window.") alias make_moving_length make_moving_sum add_help(:moving_length, "moving_length(gen, y) \ Returns the length of the values in a window over the last few inputs.") alias moving_length moving_rms # harmonicizer # (each harmonic is split into a set of harmonics via Chebyshev polynomials) # obviously very similar to ssb_bank above, but splits harmonics # individually, rather than pitch-shifting them add_help(:harmonicizer, "harmonicizer(freq, coeffs, pairs, order=40, bw=50.0, \ beg=0, dur=false, snd=false, chn=false, edpos=false) \ Splits out each harmonic and replaces it with the spectrum given in coeffs.") def harmonicizer(freq, coeffs, pairs, order = 40, bw = 50.0, beg = 0, dur = false, snd = false, chn = false, edpos = false) bands = make_array(pairs) pcoeffs = partials2polynomial(coeffs) avgs = make_array(pairs) peaks = make_array(pairs) flt = make_filter(2, vct(1, -1), vct(0, -0.9)) old_mx = maxamp new_mx = 0.0 ctr = 40 1.upto(pairs) do |i| aff = i * freq bwf = bw * (1.0 + i / (2 * pairs)) peaks[i - 1] = make_moving_max(128) avgs[i - 1] = make_moving_average(128) bands[i - 1] = make_bandpass(hz_to_2pi(aff - bwf), hz_to_2pi(aff + bwf), order) end as_one_edit_rb do map_channel_rb(beg, dur, snd, chn, edpos) do |y| sum = 0.0 bands.zip(peaks, avgs) do |bs, ps, as| sig = bandpass(bs, y) mx = moving_max(ps, sig) amp = moving_average(as, mx > 0.0 ? [100.0, 1.0 / mx].min : 0.0) if amp > 0.0 sum = sum + mx * polynomial(pcoeffs, amp * sig) end end val = filter(flt, sum) new_mx = [new_mx, val.abs].max if ctr.zero? val else ctr -= 1 0.0 end end if new_mx > 0.0 scale_channel(old_mx / new_mx, beg, dur, snd, chn) end end end # linear sampling rate conversion add_help(:linear_src_channel, "linear_src_channel(sr, snd=false, chn=false) \ Performs sampling rate conversion using linear interpolation.") def linear_src_channel(srinc, snd = false, chn = false) rd = make_sampler(0, snd, chn) last = rd.call nxt = rd.call intrp = 0.0 tempfile = with_sound(:clm, true, :output, snd_tempnam, :srate, srate(snd)) do samp = 0 until sampler_at_end?(rd) if (pos = intrp) >= 1.0 pos.floor.times do |i| last, nxt = nxt, rd.call end pos -= pos.floor end intrp = pos + pos.floor out_any(samp, last + pos * (nxt - last), 0, $output) samp += 1 end end.output len = mus_sound_framples(tempfile) set_samples(0, len - 1, tempfile, snd, chn, true, "%s(%s", get_func_name, srinc, 0, false, true) # first true=truncate to new length, false=at current edpos, # true=auto delete temp file end # Mathews/Smith High-Q filter as described in # http://ccrma.stanford.edu/~jos/smac03maxjos/ class Mfilter < Musgen def initialize(decay, freq) super() @decay = decay.to_f @frequency = freq.to_f @eps = 2.0 * sin((PI * freq) / mus_srate()) @xn = @yn = 0.0 end attr_accessor :decay, :eps def inspect format("%s.new(%0.3f, %0.3f)", self.class, @decay, @frequency) end def to_s format("#<%s decay: %0.3f, frequency: %0.3f>", self.class, @decay, @frequency) end def mfilter(x_input = 0.0, y_input = 0.0) @xn = x_input + @decay * (@xn - @eps * @yn) @yn = y_input + @decay * (@eps * @xn + @yn) end end def make_mfilter(*args) Mfilter.new(get_args(args, :decay, 0.99), get_args(args, :frequency, 1000.0)) end def mfilter(m, x_input = 0.0, y_input = 0.0) m.mfilter(x_input, y_input) end =begin with_sound(:clm, true) do rd = make_sampler(0, "now.snd") m = make_mfilter 10000.times do |i| outa(i, mfilter(m, 0.1 * rd.call), $output) end end =end # # sweep center freq: =begin with_sound(:clm, true) do rd = make_sampler(0, "oboe.snd") m = make_mfilter(:decay, 0.99, :frequency, 1000) e = make_env([0, 100, 1, 2000], :length, 10000) 10000.times do |i| outa(i, mfilter(m, 0.1 * rd.call), $output) m.eps = 2.0 * sind((PI * env(e)) / mus_srate()) end end =end # # harmonics: =begin with_sound(:clm, true, :statistics, true) do noi = make_rand(10000) filters = make_array(9) do make_mfilter(:decay, 0.999, :frequency, 400.0 * (i + 1)) end 10000.times do |i| sum = 0.0 input = 0.01 * rand(noi) filters.each do |f| sum = sum + (1.0 / (j + 1)) * mfilter(f, input) end outa(i, sum $output) end end =end # # spectrum displayed in various frequency scales # class Display_bark_fft # click in lisp-graph to change the tick placement choice def initialize @bark_fft_size = 0 @bark_tick_function = 0 end attr_reader :bark_tick_function def display_bark_fft(snd, chn) ls = left_sample(snd, chn) rs = right_sample(snd, chn) fftlen = (2 ** (log((rs - ls) + 1.0) / log(2.0)).ceil.to_i).to_i if fftlen > 0 data = channel2vct(ls, fftlen, snd, chn) normalized = (transform_normalization(snd, chn) != Dont_normalize) linear = true if vct?(data) fft = snd_spectrum(data, fft_window(snd, chn), fftlen, linear, fft_window_beta(snd, chn), false, normalized) if vct?(fft) sr = srate(snd) mx = fft.peak data_len = fft.length # bark settings bark_low = bark(20.0).floor bark_high = bark(0.5 * sr).ceil bark_frqscl = data_len / (bark_high - bark_low) bark_data = Vct.new(data_len) # mel settings mel_low = mel(20.0).floor mel_high = mel(0.5 * sr).ceil mel_frqscl = data_len / (mel_high - mel_low) mel_data = Vct.new(data_len) # erb settings erb_low = erb(20.0).floor erb_high = erb(0.5 * sr).ceil erb_frqscl = data_len / (erb_high - erb_low) erb_data = Vct.new(data_len) @bark_fft_size = fftlen fftlenf = fftlen.to_f fft.each_with_index do |val, i| frq = sr * (i / fftlenf) bark_bin = (bark_frqscl * (bark(frq) - bark_low)).round mel_bin = (mel_frqscl * (mel(frq) - mel_low)).round erb_bin = (erb_frqscl * (erb(frq) - erb_low)).round if bark_bin.between?(0, data_len - 1) bark_data[bark_bin] += val end if mel_bin.between?(0, data_len - 1) mel_data[mel_bin] += val end if erb_bin.between?(0, data_len - 1) erb_data[erb_bin] += val end end if normalized bmx = bark_data.peak mmx = mel_data.peak emx = erb_data.peak if (mx - bmx).abs > 0.01 bark_data.scale!(mx / bmx) end if (mx - mmx).abs > 0.01 mel_data.scale!(mx / mmx) end if (mx - emx).abs > 0.01 erb_data.scale!(mx / emx) end end graph([bark_data, mel_data, erb_data], "ignored", 20.0, 0.5 * sr, 0.0, (normalized ? 1.0 : data_len * y_zoom_slider(snd, chn)), snd, chn, false, Show_bare_x_axis) end end end false end def mark_bark_labels(snd, chn) # at this point the x axis has no markings, but there is room # for labels and ticks old_foreground_color = foreground_color(snd, chn, Copy_context) # assume at start the foreground color is correct axinfo = axis_info(snd, chn, Lisp_graph) axis_x0 = axinfo[10] axis_x1 = axinfo[12] axis_y0 = axinfo[13] axis_y1 = axinfo[11] label_height = 15 char_width = 8 sr2 = 0.5 * srate(snd) minor_tick_len = 6 major_tick_len = 12 tick_y0 = axis_y1 minor_y0 = axis_y1 + minor_tick_len major_y0 = axis_y1 + major_tick_len bark_label_font = snd_font(3) bark_numbers_font = snd_font(2) label_pos = (axis_x0 + 0.45 * (axis_x1 - axis_x0)).to_i cr = channel_widgets(snd, chn)[17] scale_position = lambda do |scale, f| b20 = scale.call(20.0) (axis_x0 + ((axis_x1 - axis_x0) * (scale.call(f) - b20)) / (scale.call(sr2) - b20)).round end bark_position = lambda do |f| scale_position.call(method(:bark).to_proc, f) end mel_position = lambda do |f| scale_position.call(method(:mel).to_proc, f) end erb_position = lambda do |f| scale_position.call(method(:erb).to_proc, f) end draw_bark_ticks = lambda do |bark_function| if bark_numbers_font set_current_font(bark_numbers_font, snd, chn, Copy_context) end draw_line(axis_x0, tick_y0, axis_x0, major_y0, snd, chn, Copy_context, cr) i1000 = scale_position.call(bark_function, 1000.0) i10000 = scale_position.call(bark_function, 10000.0) draw_line(i1000, tick_y0, i1000, major_y0, snd, chn, Copy_context, cr) draw_line(i10000, tick_y0, i10000, major_y0, snd, chn, Copy_context, cr) draw_string("20", axis_x0, major_y0, snd, chn, Copy_context, cr) draw_string("1000", i1000 - 3 * 4, major_y0, snd, chn, Copy_context, cr) draw_string("10000", i10000 - 6 * 4, major_y0, snd, chn, Copy_context, cr) draw_string("fft size: #{@bark_fft_size}", axis_x0 + 10, axis_y0, snd, chn, Copy_context, cr) 100.step(1000, 100) do |i| i100 = scale_position.call(bark_function, i) draw_line(i100, tick_y0, i100, minor_y0, snd, chn, Copy_context, cr) end 2000.step(10000, 1000) do |i| i1000 = scale_position.call(bark_function, i) draw_line(i1000, tick_y0, i1000, minor_y0, snd, chn, Copy_context, cr) end end # bark label/ticks if self.bark_tick_function.zero? draw_bark_ticks.call(bark_position) end if bark_label_font set_current_font(bark_label_font, snd, chn, Copy_context) end draw_string("bark,", label_pos, axis_y1 + label_height, snd, chn, Copy_context, cr) # mel label/ticks set_foreground_color(snd_color(2), snd, chn, Copy_context) if self.bark_tick_function == 1 draw_bark_ticks.call(mel_position) end if bark_label_font set_current_font(bark_label_font, snd, chn, Copy_context) end draw_string("mel,", char_width * 6 + label_pos, axis_y1 + label_height, snd, chn, Copy_context, cr) # erb label/ticks set_foreground_color(snd_color(4), snd, chn, Copy_context) if self.bark_tick_function == 2 draw_bark_ticks.call(erb_position) end if bark_label_font set_current_font(bark_label_font, snd, chn, Copy_context) end draw_string("erb", char_width * (6 + 5) + label_pos, axis_y1 + label_height, snd, chn, Copy_context, cr) set_foreground_color(old_foreground_color, snd, chn, Copy_context) end # mouse click = move to next scale's ticks def choose_bark_ticks(snd, chn, button, state, x, y, axis) if axis == Lisp_graph @bark_tick_function += 1 if @bark_tick_function > 2 @bark_tick_function = 0 end update_lisp_graph(snd, chn) end end private def bark(f) f2 = f / 7500.0 13.5 * atan(0.00076 * f) + (3.5 * atan(f2 * f2)) end def mel(f) 1127.0 * log(1.0 + f / 700.0) end def erb(f) 43.0 + 11.17 * log((f + 312.0) / (f + 14675.0)) end end # user's view of display-bark-fft function def display_bark_fft(off = false) if off.kind_of?(FalseClass) db = Display_bark_fft.new $lisp_graph_hook.add_hook!("display-bark-fft") do |snd, chn| db.display_bark_fft(snd, chn) end $after_lisp_graph_hook.add_hook!("make-bark-label") do |snd, chn| db.mark_bark_labels(snd, chn) end $mouse_click_hook.add_hook!("choose-bark-ticks") do |s, c, b, st, x, y, a| db.choose_bark_ticks(s, c, b, st, x, y, a) end Snd.sounds.each do |snd| channels(snd).times do |chn| update_lisp_graph(snd, chn) end end else $lisp_graph_hook.remove_hook!("display-bark-fft") $after_lisp_graph_hook.remove_hook!("make-bark-label") $mouse_click_hook.remove_hook!("choose-bark-ticks") Snd.sounds.each do |snd| channels(snd).times do |chn| set_lisp_graph?(false, snd, chn) end end end off end def undisplay_bark_fft display_bark_fft(true) end end include Dsp # dsp.rb ends here