diff options
Diffstat (limited to 'filters.lua')
| -rw-r--r-- | filters.lua | 361 |
1 files changed, 361 insertions, 0 deletions
diff --git a/filters.lua b/filters.lua new file mode 100644 index 0000000..e9b2900 --- /dev/null +++ b/filters.lua @@ -0,0 +1,361 @@ +--[==[ +-- +-- Non-working FIR filters (FIXME) +-- + +-- Normalized Sinc function +local function Sinc(x) + return x == 0 and 1 or + math.sin(2*math.pi*x)/(2*math.pi*x) +end + +local function Hamming(n, window) + local alpha = 0.54 + return alpha - (1-alpha)*math.cos((2*math.pi*n)/(window-1)) +end +local function Blackman(n, window) + local alpha = 0.16 + return (1-alpha)/2 - + 0.5*math.cos((2*math.pi*n)/(window-1)) + + alpha*0.5*math.cos((4*math.pi*n)/(window-1)) +end + +FIRStream = DeriveClass(Stream) + +function FIRStream:ctor(stream, freq_stream) + self.stream = tostream(stream) + self.freq_stream = tostream(freq_stream) +end + +function FIRStream:gtick() + local window = {} + + -- window size (max. 1024 samples) + -- this is the max. latency introduced by the filter + -- since the window must be filled before we can generate + -- (filtered) samples + local window_size = math.min(1024, self.stream:len()) + local window_p = window_size-1 + local accu = 0 + + local blackman = {} + for i = 1, window_size do blackman[i] = Blackman(i-1, window_size) end + + local tick = self.stream:gtick() + local freq_tick = self.freq_stream:gtick() + + return function() + -- fill buffer (initial) + while #window < window_size-1 do + table.insert(window, tick()) + end + + window[window_p+1] = tick() + window_p = (window_p + 1) % window_size + + local period = freq_tick()/samplerate + + local sample = 0 + local i = window_p + repeat + -- FIXME + sample = sample + window[(i % window_size)+1] * + Sinc((i-window_p - window_size/2)/period) * + blackman[i-window_p+1] + i = i + 1 + until (i % window_size) == window_p + + return sample + end +end + +function FIRStream:len() + return self.stream:len() +end +]==] + +-- +-- General-purpose IIR filters: +-- These are direct translations of ChucK's LPF, HPF, BPF and BRF +-- ugens which are in turn adapted from SuperCollider 3. +-- + +-- De-denormalize function adapted from ChucK. +-- Not quite sure why this is needed - properly to make the +-- IIR filters numerically more stable. +local function ddn(f) + return f >= 0 and (f > 1e-15 and f < 1e15 and f or 0) or + (f < -1e-15 and f > -1e15 and f or 0) +end + +LPFStream = DeriveClass(MuxableStream) + +function LPFStream:muxableCtor(stream, freq) + self.stream = stream + self.freq_stream = freq +end + +function LPFStream:gtick() + local a0, b1, b2 + local y1, y2 = 0, 0 + + -- some cached constants + local radians_per_sample = (2*math.pi)/samplerate + local sqrt2 = math.sqrt(2) + + -- some cached math table lookups + local tan = math.tan + + local tick = self.stream:gtick() + local freq_tick = self.freq_stream:gtick() + local cur_freq = nil + + return function() + local sample = tick() + local freq = freq_tick() + + if sample == nil or freq == nil then + -- don't filter if we run out of frequency samples + return sample + elseif freq ~= cur_freq then + -- calculate filter coefficients + -- avoid recalculation for constant frequencies + cur_freq = freq + + local pfreq = cur_freq * radians_per_sample * 0.5 + + local C = 1/tan(pfreq) + local C2 = C*C + local sqrt2C = C * sqrt2 + + a0 = 1/(1 + sqrt2C + C2) + b1 = -2.0 * (1.0 - C2) * a0 + b2 = -(1.0 - sqrt2C + C2) * a0 + end + + local y0 = sample + b1*y1 + b2*y2 + local result = a0 * (y0 + 2*y1 + y2) + + y2 = ddn(y1) + y1 = ddn(y0) + + return result + end +end + +function LPFStream:len() + return self.stream:len() +end + +function Stream:LPF(freq) + return LPFStream:new(self, freq) +end + +HPFStream = DeriveClass(MuxableStream) + +function HPFStream:muxableCtor(stream, freq) + self.stream = stream + self.freq_stream = freq +end + +function HPFStream:gtick() + local a0, b1, b2 + local y1, y2 = 0, 0 + + -- some cached constants + local radians_per_sample = (2*math.pi)/samplerate + local sqrt2 = math.sqrt(2) + + -- some cached math table lookups + local tan = math.tan + + local tick = self.stream:gtick() + local freq_tick = self.freq_stream:gtick() + local cur_freq = nil + + -- NOTE: Very similar to LPFStream.gtick() + -- Can we factor out the similarity without sacrificing + -- too much performance? + return function() + local sample = tick() + local freq = freq_tick() + + if sample == nil or freq == nil then + -- don't filter if we run out of frequency samples + return sample + elseif freq ~= cur_freq then + -- calculate filter coefficients + -- avoid recalculation for constant frequencies + cur_freq = freq + + local pfreq = cur_freq * radians_per_sample * 0.5 + + local C = tan(pfreq) + local C2 = C*C + local sqrt2C = C * sqrt2 + + a0 = 1/(1 + sqrt2C + C2) + b1 = 2.0 * (1.0 - C2) * a0 + b2 = -(1.0 - sqrt2C + C2) * a0 + end + + local y0 = sample + b1*y1 + b2*y2 + local result = a0 * (y0 - 2*y1 + y2) + + y2 = ddn(y1) + y1 = ddn(y0) + + return result + end +end + +function HPFStream:len() + return self.stream:len() +end + +function Stream:HPF(freq) + return HPFStream:new(self, freq) +end + +-- NOTE: The quality factor, indirectly proportional +-- to the passband width +BPFStream = DeriveClass(MuxableStream) + +function BPFStream:muxableCtor(stream, freq, quality) + self.stream = stream + self.freq_stream = freq + self.quality_stream = quality +end + +function BPFStream:gtick() + local a0, b1, b2 + local y1, y2 = 0, 0 + + -- some cached constants + local radians_per_sample = (2*math.pi)/samplerate + local sqrt2 = math.sqrt(2) + + -- some cached math table lookups + local tan = math.tan + local cos = math.cos + + local tick = self.stream:gtick() + local freq_tick = self.freq_stream:gtick() + local quality_tick = self.quality_stream:gtick() + local cur_freq, cur_quality + + return function() + local sample = tick() + local freq = freq_tick() + local quality = quality_tick() + + if sample == nil or freq == nil or quality == nil then + -- don't filter if we run out of frequency samples + return sample + elseif freq ~= cur_freq or quality ~= cur_quality then + -- calculate filter coefficients + -- avoid recalculation for constant frequencies + -- and quality factors + cur_freq = freq + cur_quality = quality + + local pfreq = cur_freq * radians_per_sample + local pbw = 1 / cur_quality*pfreq*0.5 + + local C = 1/tan(pbw) + local D = 2*cos(pfreq); + + a0 = 1/(1 + C) + b1 = C*D*a0 + b2 = (1 - C)*a0 + end + + local y0 = sample + b1*y1 + b2*y2 + local result = a0 * (y0 - y2) + + y2 = ddn(y1) + y1 = ddn(y0) + + return result + end +end + +function BPFStream:len() + return self.stream:len() +end + +function Stream:BPF(freq, quality) + return BPFStream:new(self, freq, quality) +end + +-- NOTE: The quality factor, indirectly proportional +-- to the passband width +BRFStream = DeriveClass(MuxableStream) + +function BRFStream:muxableCtor(stream, freq, quality) + self.stream = stream + self.freq_stream = freq + self.quality_stream = quality +end + +function BRFStream:gtick() + local a0, b1, b2 + local y1, y2 = 0, 0 + + -- some cached constants + local radians_per_sample = (2*math.pi)/samplerate + local sqrt2 = math.sqrt(2) + + -- some cached math table lookups + local tan = math.tan + local cos = math.cos + + local tick = self.stream:gtick() + local freq_tick = self.freq_stream:gtick() + local quality_tick = self.quality_stream:gtick() + local cur_freq, cur_quality + + -- NOTE: Very similar to BPFStream.gtick() + return function() + local sample = tick() + local freq = freq_tick() + local quality = quality_tick() + + if sample == nil or freq == nil or quality == nil then + -- don't filter if we run out of frequency samples + return sample + elseif freq ~= cur_freq then + -- calculate filter coefficients + -- avoid recalculation for constant frequencies + -- and quality factors + cur_freq = freq + cur_quality = quality + + local pfreq = cur_freq * radians_per_sample + local pbw = 1 / cur_quality*pfreq*0.5 + + local C = tan(pbw) + local D = 2*cos(pfreq); + + a0 = 1/(1 + C) + b1 = -D*a0 + b2 = (1 - C)*a0 + end + + local y0 = sample - b1*y1 - b2*y2 + local result = a0 * (y0 + y2) + b1*y1 + + y2 = ddn(y1) + y1 = ddn(y0) + + return result + end +end + +function BRFStream:len() + return self.stream:len() +end + +function Stream:BRF(freq, quality) + return BRFStream:new(self, freq, quality) +end |