diff --git a/OSBindings/Mac/Clock Signal.xcodeproj/project.pbxproj b/OSBindings/Mac/Clock Signal.xcodeproj/project.pbxproj index 9ba0da4f8..6ea4d6d7a 100644 --- a/OSBindings/Mac/Clock Signal.xcodeproj/project.pbxproj +++ b/OSBindings/Mac/Clock Signal.xcodeproj/project.pbxproj @@ -308,6 +308,7 @@ 4BBF99161C8FBA6F0075DAFB /* CRTRunBuilder.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 4BBF990C1C8FBA6F0075DAFB /* CRTRunBuilder.cpp */; }; 4BBF99171C8FBA6F0075DAFB /* Shader.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 4BBF99101C8FBA6F0075DAFB /* Shader.cpp */; }; 4BBF99181C8FBA6F0075DAFB /* TextureTarget.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 4BBF99121C8FBA6F0075DAFB /* TextureTarget.cpp */; }; + 4BC76E691C98E31700E6EF73 /* FIRFilter.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 4BC76E671C98E31700E6EF73 /* FIRFilter.cpp */; }; 4BCB70B41C947DDC005B1712 /* plus1.rom in Resources */ = {isa = PBXBuildFile; fileRef = 4BCB70B31C947DDC005B1712 /* plus1.rom */; }; 4BE5F85E1C3E1C2500C43F01 /* basic.rom in Resources */ = {isa = PBXBuildFile; fileRef = 4BE5F85C1C3E1C2500C43F01 /* basic.rom */; }; 4BE5F85F1C3E1C2500C43F01 /* os.rom in Resources */ = {isa = PBXBuildFile; fileRef = 4BE5F85D1C3E1C2500C43F01 /* os.rom */; }; @@ -666,6 +667,8 @@ 4BBF99121C8FBA6F0075DAFB /* TextureTarget.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = TextureTarget.cpp; sourceTree = ""; }; 4BBF99131C8FBA6F0075DAFB /* TextureTarget.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = TextureTarget.hpp; sourceTree = ""; }; 4BBF99191C8FC2750075DAFB /* CRTTypes.hpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.h; path = CRTTypes.hpp; sourceTree = ""; }; + 4BC76E671C98E31700E6EF73 /* FIRFilter.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = FIRFilter.cpp; sourceTree = ""; }; + 4BC76E681C98E31700E6EF73 /* FIRFilter.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = FIRFilter.hpp; sourceTree = ""; }; 4BCB70B31C947DDC005B1712 /* plus1.rom */ = {isa = PBXFileReference; lastKnownFileType = file; path = plus1.rom; sourceTree = ""; }; 4BE5F85C1C3E1C2500C43F01 /* basic.rom */ = {isa = PBXFileReference; lastKnownFileType = file; path = basic.rom; sourceTree = ""; }; 4BE5F85D1C3E1C2500C43F01 /* os.rom */ = {isa = PBXFileReference; lastKnownFileType = file; path = os.rom; sourceTree = ""; }; @@ -733,6 +736,8 @@ 4B2409591C45DF85004DA684 /* SignalProcessing */ = { isa = PBXGroup; children = ( + 4BC76E671C98E31700E6EF73 /* FIRFilter.cpp */, + 4BC76E681C98E31700E6EF73 /* FIRFilter.hpp */, 4B24095A1C45DF85004DA684 /* Stepper.hpp */, ); name = SignalProcessing; @@ -1646,6 +1651,7 @@ 4B55CE581C3B7D360093A61B /* Atari2600Document.swift in Sources */, 4B0EBFB81C487F2F00A11F35 /* AudioQueue.m in Sources */, 4BBF99181C8FBA6F0075DAFB /* TextureTarget.cpp in Sources */, + 4BC76E691C98E31700E6EF73 /* FIRFilter.cpp in Sources */, 4B55CE5F1C3B7D960093A61B /* MachineDocument.swift in Sources */, 4B69FB441C4D941400B5F0AA /* TapeUEF.cpp in Sources */, 4BBF99141C8FBA6F0075DAFB /* CRTInputBufferBuilder.cpp in Sources */, diff --git a/SignalProcessing/FIRFilter.cpp b/SignalProcessing/FIRFilter.cpp new file mode 100644 index 000000000..2668f8a55 --- /dev/null +++ b/SignalProcessing/FIRFilter.cpp @@ -0,0 +1,148 @@ +// +// LinearFilter.c +// Clock Signal +// +// Created by Thomas Harte on 01/10/2011. +// Copyright 2011 Thomas Harte. All rights reserved. +// + +#include "FIRFilter.hpp" +#include + +using namespace SignalProcessing; + +/* + + A Kaiser-Bessel filter is a real time window filter. It looks at the last n samples + of an incoming data source and computes a filtered value, which is the value you'd + get after applying the specified filter, at the centre of the sampling window. + + Hence, if you request a 37 tap filter then filtering introduces a latency of 18 + samples. Suppose you're receiving input at 44,100Hz and using 4097 taps, then you'll + introduce a latency of 2048 samples, which is about 46ms. + + There's a correlation between the number of taps and the quality of the filtering. + More samples = better filtering, at the cost of greater latency. Internally, applying + the filter involves calculating a weighted sum of previous values, so increasing the + number of taps is quite cheap in processing terms. + + Original source for this filter: + + "DIGITAL SIGNAL PROCESSING, II", IEEE Press, pages 123–126. +*/ + + +// our little fixed point scheme +#define kCSKaiserBesselFilterFixedMultiplier 32767.0f +#define kCSKaiserBesselFilterFixedShift 15 + +/* ino evaluates the 0th order Bessel function at a */ +static float csfilter_ino(float a) +{ + float d = 0.0f; + float ds = 1.0f; + float s = 1.0f; + + do + { + d += 2.0f; + ds *= (a * a) / (d * d); + s += ds; + } + while(ds > s*1e-6f); + + return s; +} + +static void csfilter_setIdealisedFilterResponse(short *filterCoefficients, float *A, float attenuation, unsigned int numberOfTaps) +{ + /* calculate alpha, which is the Kaiser-Bessel window shape factor */ + float a; // to take the place of alpha in the normal derivation + + if(attenuation < 21.0f) + a = 0.0f; + else + { + if(attenuation > 50.0f) + a = 0.1102f * (attenuation - 8.7f); + else + a = 0.5842f * powf(attenuation - 21.0f, 0.4f) + 0.7886f * (attenuation - 21.0f); + } + + float *filterCoefficientsFloat = (float *)malloc(sizeof(float) * numberOfTaps); + + /* work out the right hand side of the filter coefficients */ + unsigned int Np = (numberOfTaps - 1) / 2; + float I0 = csfilter_ino(a); + float NpSquared = (float)(Np * Np); + for(unsigned int i = 0; i <= Np; i++) + { + filterCoefficientsFloat[Np + i] = + A[i] * + csfilter_ino(a * sqrtf(1.0f - ((float)(i * i) / NpSquared) )) / + I0; + } + + /* coefficients are symmetrical, so copy from right hand side to left side */ + for(unsigned int i = 0; i < Np; i++) + { + filterCoefficientsFloat[i] = filterCoefficientsFloat[numberOfTaps - 1 - i]; + } + + /* scale back up so that we retain 100% of input volume */ + float coefficientTotal = 0.0f; + for(unsigned int i = 0; i < numberOfTaps; i++) + { + coefficientTotal += filterCoefficientsFloat[i]; + } + + /* we'll also need integer versions, potentially */ + float coefficientMultiplier = 1.0f / coefficientTotal; + for(unsigned int i = 0; i < numberOfTaps; i++) + { + filterCoefficients[i] = (short)(filterCoefficientsFloat[i] * kCSKaiserBesselFilterFixedMultiplier * coefficientMultiplier); + } + + free(filterCoefficientsFloat); +} + +FIRFilter::FIRFilter(unsigned int number_of_taps, unsigned int input_sample_rate, float low_frequency, float high_frequency, float attenuation) +{ + // we must be asked to filter based on an odd number of + // taps, and at least three + if(number_of_taps < 3) number_of_taps = 3; + if(attenuation < 21.0f) attenuation = 21.0f; + + // ensure we have an odd number of taps + number_of_taps |= 1; + + // store instance variables + number_of_taps_ = number_of_taps; + filter_coefficients_ = new short[number_of_taps_]; + + /* calculate idealised filter response */ + unsigned int Np = (number_of_taps - 1) / 2; + float twoOverSampleRate = 2.0f / (float)input_sample_rate; + + float *A = new float[Np+1]; + A[0] = 2.0f * (high_frequency - low_frequency) / (float)input_sample_rate; + for(unsigned int i = 1; i <= Np; i++) + { + float iPi = (float)i * (float)M_PI; + A[i] = + ( + sinf(twoOverSampleRate * iPi * high_frequency) - + sinf(twoOverSampleRate * iPi * low_frequency) + ) / iPi; + } + + csfilter_setIdealisedFilterResponse(filter_coefficients_, A, attenuation, number_of_taps_); + + /* clean up */ + delete[] A; +} + +FIRFilter::~FIRFilter() +{ + delete[] filter_coefficients_; +} diff --git a/SignalProcessing/FIRFilter.hpp b/SignalProcessing/FIRFilter.hpp new file mode 100644 index 000000000..27319d055 --- /dev/null +++ b/SignalProcessing/FIRFilter.hpp @@ -0,0 +1,83 @@ +// +// LinearFilter.h +// Clock Signal +// +// Created by Thomas Harte on 01/10/2011. +// Copyright 2011 Thomas Harte. All rights reserved. +// + +#ifndef FIRFilter_hpp +#define FIRFilter_hpp + +/* + + The FIR filter takes a 1d PCM signal with + a given sample rate and filters it according + to a specified filter (band pass only at + present, more to come if required). The number + of taps (ie, samples considered simultaneously + to make an output sample) is configurable; + smaller numbers permit a filter that operates + more quickly and with less lag but less + effectively. + + FIR filters are window functions; expected use is + to point sample an input that has been subject to + a filter. + +*/ + +#ifdef __APPLE__ +#include +#endif + +namespace SignalProcessing { + +class FIRFilter { + public: + /*! + Creates an instance of @c FIRFilter. + + @param number_of_taps The size of window for input data. + @param input_sample_rate The sampling rate of the input signal. + @param low_frequency The lowest frequency of signal to retain in the output. + @param high_frequency The highest frequency of signal to retain in the output. + @param attenuation The attenuation of the output. + */ + FIRFilter(unsigned int number_of_taps, unsigned int input_sample_rate, float low_frequency, float high_frequency, float attenuation); + + ~FIRFilter(); + + /*! A suggested default attenuation value. */ + const float DefaultAttenuation = 60.0f; + + /*! + Applies the filter to one batch of input samples, returning the net result. + + @param src The source buffer to apply the filter to. + @returns The result of applying the filter. + */ + inline short apply(const short *src) + { + #ifdef __APPLE__ + short result; + vDSP_dotpr_s1_15(filter_coefficients_, 1, src, 1, &result, number_of_taps_); + return result; + #else + int outputValue = 0; + for(unsigned int c = 0; c < number_of_taps_; c++) + { + outputValue += filter_coefficients_[c] * src[c]; + } + return (short)(outputValue >> kCSKaiserBesselFilterFixedShift); + #endif + } + + private: + short *filter_coefficients_; + unsigned int number_of_taps_; +}; + +} + +#endif