001/**
002 * Copyright (c) 2011, The University of Southampton and the individual contributors.
003 * All rights reserved.
004 *
005 * Redistribution and use in source and binary forms, with or without modification,
006 * are permitted provided that the following conditions are met:
007 *
008 *   *  Redistributions of source code must retain the above copyright notice,
009 *      this list of conditions and the following disclaimer.
010 *
011 *   *  Redistributions in binary form must reproduce the above copyright notice,
012 *      this list of conditions and the following disclaimer in the documentation
013 *      and/or other materials provided with the distribution.
014 *
015 *   *  Neither the name of the University of Southampton nor the names of its
016 *      contributors may be used to endorse or promote products derived from this
017 *      software without specific prior written permission.
018 *
019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
029 */
030/**
031 * 
032 */
033package org.openimaj.audio.filters;
034
035import org.openimaj.audio.AudioStream;
036import org.openimaj.audio.SampleChunk;
037import org.openimaj.audio.processor.AudioProcessor;
038import org.openimaj.audio.samples.SampleBuffer;
039
040/**
041 *      A class that encapsulates a bunch of different EQ-based filter algorithms
042 *      that use a standard bi-quad filter (4th order Linkwitz-Riley filter).
043 *
044 *      @see "http://musicdsp.org/showArchiveComment.php?ArchiveID=266"
045 *      @author David Dupplaw (dpd@ecs.soton.ac.uk)
046 *  @created 19 Jul 2012
047 *      @version $Author$, $Revision$, $Date$
048 */
049public class EQFilter extends AudioProcessor
050{
051        /**
052         *      A class that represents the bi-quad coefficients for a filter.
053         *
054         *      @author David Dupplaw (dpd@ecs.soton.ac.uk)
055         *  @created 20 Jul 2012
056         *      @version $Author$, $Revision$, $Date$
057         */
058        public static class EQVars
059        {
060                double a0, a1, a2, a3, a4;
061                double b0, b1, b2, b3, b4;
062        }
063        
064        /**
065         *      An enumerator for various audio filters.
066         *
067         *      @author David Dupplaw (dpd@ecs.soton.ac.uk)
068         *  @created 19 Jul 2012
069         *      @version $Author$, $Revision$, $Date$
070         */
071        public static enum EQType
072        {
073                /**
074                 *      Low pass filter 
075                 */
076                LPF
077                {
078                        @Override
079                        public EQVars getCoefficients( final double frequency, final double sampleRate )
080                        {
081                                final double wc  = 2 * Math.PI * frequency;
082                                final double wc2 = wc * wc;
083                                final double wc3 = wc2 * wc;
084                                final double wc4 = wc2 * wc2;
085                                final double k = wc/Math.tan(Math.PI*frequency/sampleRate);
086                                final double k2 = k * k;
087                                final double k3 = k2 * k;
088                                final double k4 = k2 * k2;
089                                final double sqrt2 = Math.sqrt(2);
090                                final double sq_tmp1 = sqrt2 * wc3 * k;
091                                final double sq_tmp2 = sqrt2 * wc * k3;
092                                final double a_tmp = 4 * wc2 * k2 + 2*sq_tmp1 + k4 + 2*sq_tmp2 + wc4;
093
094                                final EQVars v = new EQVars();
095                                
096                                v.b1 = (4*(wc4+sq_tmp1-k4-sq_tmp2)) / a_tmp;
097                                v.b2 = (6*wc4-8*wc2*k2+6*k4)/a_tmp;
098                                v.b3 = (4*(wc4-sq_tmp1+sq_tmp2-k4))/a_tmp;
099                                v.b4 = (k4-2*sq_tmp1+wc4-2*sq_tmp2+4*wc2*k2)/a_tmp;
100
101                                //================================================
102                                // low-pass
103                                //================================================
104                                v.a0 = wc4/a_tmp;
105                                v.a1 = 4*wc4/a_tmp;
106                                v.a2 = 6*wc4/a_tmp;
107                                v.a3 = v.a1;
108                                v.a4 = v.a0;
109                                
110                                return v;
111                        }
112                },
113                
114                /**
115                 *      High pass filter 
116                 */
117                HPF
118                {
119                        @Override
120                        public EQVars getCoefficients( final double frequency, final double sampleRate )
121                        {
122                                final double wc  = 2 * Math.PI * frequency;
123                                final double wc2 = wc * wc;
124                                final double wc3 = wc2 * wc;
125                                final double wc4 = wc2 * wc2;
126                                final double k = wc/Math.tan(Math.PI*frequency/sampleRate);
127                                final double k2 = k * k;
128                                final double k3 = k2 * k;
129                                final double k4 = k2 * k2;
130                                final double sqrt2 = Math.sqrt(2);
131                                final double sq_tmp1 = sqrt2 * wc3 * k;
132                                final double sq_tmp2 = sqrt2 * wc * k3;
133                                final double a_tmp = 4 * wc2 * k2 + 2*sq_tmp1 + k4 + 2*sq_tmp2 + wc4;
134
135                                final EQVars v = new EQVars();
136                                
137                                v.b1 = (4*(wc4+sq_tmp1-k4-sq_tmp2)) / a_tmp;
138                                v.b2 = (6*wc4-8*wc2*k2+6*k4)/a_tmp;
139                                v.b3 = (4*(wc4-sq_tmp1+sq_tmp2-k4))/a_tmp;
140                                v.b4 = (k4-2*sq_tmp1+wc4-2*sq_tmp2+4*wc2*k2)/a_tmp;
141
142                                //================================================
143                                // high-pass
144                                //================================================
145                                v.a0 = k4/a_tmp;
146                                v.a1 = -4*k4/a_tmp;
147                                v.a2 = 6*k4/a_tmp;
148                                v.a3 = v.a1;
149                                v.a4 = v.a0;
150                                
151                                return v;
152                        }
153                };
154                
155                /**
156                 *      Initialise the filter
157                 *      @param frequency The frequency of the filter
158                 *      @param sampleRate The sample rate of the samples
159                 *      @return The biquad coefficients for this filter 
160                 */
161                public abstract EQVars getCoefficients( double frequency, double sampleRate );
162        }
163        
164        /** The type of EQ process */
165        private EQType type = null;
166        
167        /** The cached biquad coefficients for this filter */
168        private EQVars vars = null;
169        
170        /** The frequency of the filter */
171        private double frequency = 0;
172
173        // These variables are used during the loop
174        // and are stored between loops to avoid clicking
175        private double xm1 = 0, xm2 = 0, xm3 = 0, xm4 = 0, 
176                                   ym1 = 0, ym2 = 0, ym3 = 0, ym4 = 0;
177        
178        /**
179         *      Default constructor for ad-hoc processing.
180         *      @param type The type of EQ
181         *      @param f The frequency of the filter 
182         */
183        public EQFilter( final EQType type, final double f )
184        {
185                this.type = type;
186                this.frequency = f;
187        }
188        
189        /**
190         *      Chainable constructor for stream processing.
191         *      @param as The audio stream to process
192         *      @param type The type of EQ
193         *      @param f The frequency of the filter 
194         */
195        public EQFilter( final AudioStream as, final EQType type, final double f )
196        {
197                super( as );
198                this.type = type;
199                this.frequency = f;
200        }
201        
202        /**
203         *      {@inheritDoc}
204         *      @see org.openimaj.audio.processor.AudioProcessor#process(org.openimaj.audio.SampleChunk)
205         */
206        @Override
207        public SampleChunk process( final SampleChunk sample ) throws Exception
208        {
209                if( this.vars == null )
210                        this.vars = this.type.getCoefficients( this.frequency, 
211                                        sample.getFormat().getSampleRateKHz() *1000d );
212
213                // Standard bi-quad processing of the sample chunk
214                // We have to process each channel independently because the function
215                // is recursive.
216                final SampleBuffer sb = sample.getSampleBuffer();
217                final int nChans = sample.getFormat().getNumChannels();
218                for( int c = 0; c < nChans; c++ )
219                {
220                        this.xm1 = this.xm2 = this.xm3 = this.xm4 = 0;
221                        this.ym1 = this.ym2 = this.ym3 = this.ym4 = 0;
222                        for( int n = c; n < sb.size(); n += nChans )
223                        {
224                                final double tempx = sb.get(n);
225                                final double tempy = this.vars.a0*tempx + this.vars.a1*this.xm1 
226                                                + this.vars.a2*this.xm2 + this.vars.a3*this.xm3 + 
227                                                this.vars.a4*this.xm4 - this.vars.b1*this.ym1 - 
228                                                this.vars.b2*this.ym2 - this.vars.b3*this.ym3 - 
229                                                this.vars.b4*this.ym4;
230                                
231                                this.xm4 = this.xm3; this.xm3 = this.xm2; this.xm2 = this.xm1;
232                                this.xm1 = tempx;
233                                this.ym4 = this.ym3; this.ym3 = this.ym2; this.ym2 = this.ym1;
234                                this.ym1 = tempy;
235                
236                                sb.set( n, (float)tempy );
237                        }
238                }
239                        
240                return sample;
241        }
242}