3rdparty/SoundTouch/AAFilter.cpp

////////////////////////////////////////////////////////////////////////////////
///
/// FIR low-pass (anti-alias) filter with filter coefficient design routine and
/// MMX optimization.
///
/// Anti-alias filter is used to prevent folding of high frequencies when
/// transposing the sample rate with interpolation.
///
/// Author        : Copyright (c) Olli Parviainen
/// Author e-mail : oparviai 'at' iki.fi
/// SoundTouch WWW: http://www.surina.net/soundtouch
///
////////////////////////////////////////////////////////////////////////////////
//
// Last changed  : $Date: 2009-01-11 13:34:24 +0200 (Sun, 11 Jan 2009) $
// File revision : $Revision: 4 $
//
// $Id: AAFilter.cpp 45 2009-01-11 11:34:24Z oparviai $
//
////////////////////////////////////////////////////////////////////////////////
//
// License :
//
//  SoundTouch audio processing library
//  Copyright (c) Olli Parviainen
//
//  This library is free software; you can redistribute it and/or
//  modify it under the terms of the GNU Lesser General Public
//  License as published by the Free Software Foundation; either
//  version 2.1 of the License, or (at your option) any later version.
//
//  This library is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//  Lesser General Public License for more details.
//
//  You should have received a copy of the GNU Lesser General Public
//  License along with this library; if not, write to the Free Software
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
////////////////////////////////////////////////////////////////////////////////

#include <memory.h>
#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include "AAFilter.h"
#include "FIRFilter.h"

using namespace soundtouch;

#define PI        3.141592655357989
#define TWOPI    (2 * PI)

/*****************************************************************************
 *
 * Implementation of the class 'AAFilter'
 *
 *****************************************************************************/

AAFilter::AAFilter(uint len)
{
    pFIR = FIRFilter::newInstance();
    cutoffFreq = 0.5;
    setLength(len);
}


AAFilter::~AAFilter()
{
    delete pFIR;
}


// Sets new anti-alias filter cut-off edge frequency, scaled to
// sampling frequency (nyquist frequency = 0.5).
// The filter will cut frequencies higher than the given frequency.
void AAFilter::setCutoffFreq(double newCutoffFreq)
{
    cutoffFreq = newCutoffFreq;
    calculateCoeffs();
}


// Sets number of FIR filter taps
void AAFilter::setLength(uint newLength)
{
    length = newLength;
    calculateCoeffs();
}


// Calculates coefficients for a low-pass FIR filter using Hamming window
void AAFilter::calculateCoeffs()
{
    uint i;
    double cntTemp, temp, tempCoeff,h, w;
    double fc2, wc;
    double scaleCoeff, sum;
    double *work;
    SAMPLETYPE *coeffs;

    assert(length >= 2);
    assert(length % 4 == 0);
    assert(cutoffFreq >= 0);
    assert(cutoffFreq <= 0.5);

    work = new double[length];
    coeffs = new SAMPLETYPE[length];

    fc2 = 2.0 * cutoffFreq;
    wc = PI * fc2;
    tempCoeff = TWOPI / (double)length;

    sum = 0;
    for (i = 0; i < length; i ++)
    {
        cntTemp = (double)i - (double)(length / 2);

        temp = cntTemp * wc;
        if (temp != 0)
        {
            h = fc2 * sin(temp) / temp;                     // sinc function
        }
        else
        {
            h = 1.0;
        }
        w = 0.54 + 0.46 * cos(tempCoeff * cntTemp);       // hamming window

        temp = w * h;
        work[i] = temp;

        // calc net sum of coefficients
        sum += temp;
    }

    // ensure the sum of coefficients is larger than zero
    assert(sum > 0);

    // ensure we've really designed a lowpass filter...
    assert(work[length/2] > 0);
    assert(work[length/2 + 1] > -1e-6);
    assert(work[length/2 - 1] > -1e-6);

    // Calculate a scaling coefficient in such a way that the result can be
    // divided by 16384
    scaleCoeff = 16384.0f / sum;

    for (i = 0; i < length; i ++)
    {
        // scale & round to nearest integer
        temp = work[i] * scaleCoeff;
        temp += (temp >= 0) ? 0.5 : -0.5;
        // ensure no overfloods
        assert(temp >= -32768 && temp <= 32767);
        coeffs[i] = (SAMPLETYPE)temp;
    }

    // Set coefficients. Use divide factor 14 => divide result by 2^14 = 16384
    pFIR->setCoefficients(coeffs, length, 14);

    delete[] work;
    delete[] coeffs;
}


// Applies the filter to the given sequence of samples.
// Note : The amount of outputted samples is by value of 'filter length'
// smaller than the amount of input samples.
uint AAFilter::evaluate(SAMPLETYPE *dest, const SAMPLETYPE *src, uint numSamples, uint numChannels) const
{
    return pFIR->evaluate(dest, src, numSamples, numChannels);
}


uint AAFilter::getLength() const
{
    return pFIR->getLength();
}
1	////////////////////////////////////////////////////////////////////////////////
2	///
3	/// FIR low-pass (anti-alias) filter with filter coefficient design routine and
4	/// MMX optimization.
5	///
6	/// Anti-alias filter is used to prevent folding of high frequencies when
7	/// transposing the sample rate with interpolation.
8	///
9	/// Author : Copyright (c) Olli Parviainen
10	/// Author e-mail : oparviai 'at' iki.fi
11	/// SoundTouch WWW: http://www.surina.net/soundtouch
12	///
13	////////////////////////////////////////////////////////////////////////////////
14	//
15	// Last changed : $Date: 2009-01-11 13:34:24 +0200 (Sun, 11 Jan 2009) $
16	// File revision : $Revision: 4 $
17	//
18	// $Id: AAFilter.cpp 45 2009-01-11 11:34:24Z oparviai $
19	//
20	////////////////////////////////////////////////////////////////////////////////
21	//
22	// License :
23	//
24	// SoundTouch audio processing library
25	// Copyright (c) Olli Parviainen
26	//
27	// This library is free software; you can redistribute it and/or
28	// modify it under the terms of the GNU Lesser General Public
29	// License as published by the Free Software Foundation; either
30	// version 2.1 of the License, or (at your option) any later version.
31	//
32	// This library is distributed in the hope that it will be useful,
33	// but WITHOUT ANY WARRANTY; without even the implied warranty of
34	// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
35	// Lesser General Public License for more details.
36	//
37	// You should have received a copy of the GNU Lesser General Public
38	// License along with this library; if not, write to the Free Software
39	// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
40	//
41	////////////////////////////////////////////////////////////////////////////////
42
43	#include <memory.h>
44	#include <assert.h>
45	#include <math.h>
46	#include <stdlib.h>
47	#include "AAFilter.h"
48	#include "FIRFilter.h"
49
50	using namespace soundtouch;
51
52	#define PI 3.141592655357989
53	#define TWOPI (2 * PI)
54
55	/*****************************************************************************
56	*
57	* Implementation of the class 'AAFilter'
58	*
59	*****************************************************************************/
60
61	AAFilter::AAFilter(uint len)
62	{
63	pFIR = FIRFilter::newInstance();
64	cutoffFreq = 0.5;
65	setLength(len);
66	}
67
68
69
70	AAFilter::~AAFilter()
71	{
72	delete pFIR;
73	}
74
75
76
77	// Sets new anti-alias filter cut-off edge frequency, scaled to
78	// sampling frequency (nyquist frequency = 0.5).
79	// The filter will cut frequencies higher than the given frequency.
80	void AAFilter::setCutoffFreq(double newCutoffFreq)
81	{
82	cutoffFreq = newCutoffFreq;
83	calculateCoeffs();
84	}
85
86
87
88	// Sets number of FIR filter taps
89	void AAFilter::setLength(uint newLength)
90	{
91	length = newLength;
92	calculateCoeffs();
93	}
94
95
96
97	// Calculates coefficients for a low-pass FIR filter using Hamming window
98	void AAFilter::calculateCoeffs()
99	{
100	uint i;
101	double cntTemp, temp, tempCoeff,h, w;
102	double fc2, wc;
103	double scaleCoeff, sum;
104	double *work;
105	SAMPLETYPE *coeffs;
106
107	assert(length >= 2);
108	assert(length % 4 == 0);
109	assert(cutoffFreq >= 0);
110	assert(cutoffFreq <= 0.5);
111
112	work = new double[length];
113	coeffs = new SAMPLETYPE[length];
114
115	fc2 = 2.0 * cutoffFreq;
116	wc = PI * fc2;
117	tempCoeff = TWOPI / (double)length;
118
119	sum = 0;
120	for (i = 0; i < length; i ++)
121	{
122	cntTemp = (double)i - (double)(length / 2);
123
124	temp = cntTemp * wc;
125	if (temp != 0)
126	{
127	h = fc2 * sin(temp) / temp; // sinc function
128	}
129	else
130	{
131	h = 1.0;
132	}
133	w = 0.54 + 0.46 * cos(tempCoeff * cntTemp); // hamming window
134
135	temp = w * h;
136	work[i] = temp;
137
138	// calc net sum of coefficients
139	sum += temp;
140	}
141
142	// ensure the sum of coefficients is larger than zero
143	assert(sum > 0);
144
145	// ensure we've really designed a lowpass filter...
146	assert(work[length/2] > 0);
147	assert(work[length/2 + 1] > -1e-6);
148	assert(work[length/2 - 1] > -1e-6);
149
150	// Calculate a scaling coefficient in such a way that the result can be
151	// divided by 16384
152	scaleCoeff = 16384.0f / sum;
153
154	for (i = 0; i < length; i ++)
155	{
156	// scale & round to nearest integer
157	temp = work[i] * scaleCoeff;
158	temp += (temp >= 0) ? 0.5 : -0.5;
159	// ensure no overfloods
160	assert(temp >= -32768 && temp <= 32767);
161	coeffs[i] = (SAMPLETYPE)temp;
162	}
163
164	// Set coefficients. Use divide factor 14 => divide result by 2^14 = 16384
165	pFIR->setCoefficients(coeffs, length, 14);
166
167	delete[] work;
168	delete[] coeffs;
169	}
170
171
172	// Applies the filter to the given sequence of samples.
173	// Note : The amount of outputted samples is by value of 'filter length'
174	// smaller than the amount of input samples.
175	uint AAFilter::evaluate(SAMPLETYPE dest, const SAMPLETYPE src, uint numSamples, uint numChannels) const
176	{
177	return pFIR->evaluate(dest, src, numSamples, numChannels);
178	}
179
180
181	uint AAFilter::getLength() const
182	{
183	return pFIR->getLength();
184	}