|
|
|
|
 |
Beat Detector Class
References : Posted by DSPMaster[at]free[dot]fr
Notes :
This class was designed for a VST plugin. Basically, it's just a 2nd order LP filter, followed by an enveloppe detector (thanks Bram), feeding a Schmitt trigger. The rising edge detector provides a 1-sample pulse each time a beat is detected. Code is self documented...
Note : The class uses a fixed comparison level, you may need to change it.
Code :
// ***** BEATDETECTOR.H *****
#ifndef BeatDetectorH
#define BeatDetectorH
class TBeatDetector
{
private:
float KBeatFilter; // Filter coefficient
float Filter1Out, Filter2Out;
float BeatRelease; // Release time coefficient
float PeakEnv; // Peak enveloppe follower
bool BeatTrigger; // Schmitt trigger output
bool PrevBeatPulse; // Rising edge memory
public:
bool BeatPulse; // Beat detector output
TBeatDetector();
~TBeatDetector();
virtual void setSampleRate(float SampleRate);
virtual void AudioProcess (float input);
};
#endif
// ***** BEATDETECTOR.CPP *****
#include "BeatDetector.h"
#include "math.h"
#define FREQ_LP_BEAT 150.0f // Low Pass filter frequency
#define T_FILTER 1.0f/(2.0f*M_PI*FREQ_LP_BEAT) // Low Pass filter time constant
#define BEAT_RTIME 0.02f // Release time of enveloppe detector in second
TBeatDetector::TBeatDetector()
// Beat detector constructor
{
Filter1Out=0.0;
Filter2Out=0.0;
PeakEnv=0.0;
BeatTrigger=false;
PrevBeatPulse=false;
setSampleRate(44100);
}
TBeatDetector::~TBeatDetector()
{
// Nothing specific to do...
}
void TBeatDetector::setSampleRate (float sampleRate)
// Compute all sample frequency related coeffs
{
KBeatFilter=1.0/(sampleRate*T_FILTER);
BeatRelease=(float)exp(-1.0f/(sampleRate*BEAT_RTIME));
}
void TBeatDetector::AudioProcess (float input)
// Process incoming signal
{
float EnvIn;
// Step 1 : 2nd order low pass filter (made of two 1st order RC filter)
Filter1Out=Filter1Out+(KBeatFilter*(input-Filter1Out));
Filter2Out=Filter2Out+(KBeatFilter*(Filter1Out-Filter2Out));
// Step 2 : peak detector
EnvIn=fabs(Filter2Out);
if (EnvIn>PeakEnv) PeakEnv=EnvIn; // Attack time = 0
else
{
PeakEnv*=BeatRelease;
PeakEnv+=(1.0f-BeatRelease)*EnvIn;
}
// Step 3 : Schmitt trigger
if (!BeatTrigger)
{
if (PeakEnv>0.3) BeatTrigger=true;
}
else
{
if (PeakEnv<0.15) BeatTrigger=false;
}
// Step 4 : rising edge detector
BeatPulse=false;
if ((BeatTrigger)&&(!PrevBeatPulse))
BeatPulse=true;
PrevBeatPulse=BeatTrigger;
}
4 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
DFT
Type : fourier transform
References : Posted by Andy Mucho
Code :
AnalyseWaveform(float *waveform, int framesize)
{
float aa[MaxPartials];
float bb[MaxPartials];
for(int i=0;i<partials;i++)
{
aa[i]=0;
bb[i]=0;
}
int hfs=framesize/2;
float pd=pi/hfs;
for (i=0;i<framesize;i++)
{
float w=waveform[i];
int im = i-hfs;
for(int h=0;h<partials;h++)
{
float th=(pd*(h+1))*im;
aa[h]+=w*cos(th);
bb[h]+=w*sin(th);
}
}
for (int h=0;h<partials;h++)
amp[h]= sqrt(aa[h]*aa[h]+bb[h]*bb[h])/hfs;
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Envelope detector
References : Posted by Bram
Notes :
Basicaly a one-pole LP filter with different coefficients for attack and release fed by the abs() of the signal. If you don't need different attack and decay settings, just use in->abs()->LP
Code :
//attack and release in milliseconds
float ga = (float) exp(-1/(SampleRate*attack));
float gr = (float) exp(-1/(SampleRate*release));
float envelope=0;
for(...)
{
//get your data into 'input'
EnvIn = abs(input);
if(envelope < EnvIn)
{
envelope *= ga;
envelope += (1-ga)*EnvIn;
}
else
{
envelope *= gr;
envelope += (1-gr)*EnvIn;
}
//envelope now contains.........the envelope ;)
}
4 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Envelope follower with different attack and release
References : Posted by Bram
Notes :
xxxx_in_ms is xxxx in milliseconds ;-)
Code :
init::
attack_coef = exp(log(0.01)/( attack_in_ms * samplerate * 0.001));
release_coef = exp(log(0.01)/( release_in_ms * samplerate * 0.001));
envelope = 0.0;
loop::
tmp = fabs(in);
if(tmp > envelope)
envelope = attack_coef * (envelope - tmp) + tmp;
else
envelope = release_coef * (envelope - tmp) + tmp;
4 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Fast in-place Walsh-Hadamard Transform
Type : wavelet transform
References : Posted by Timo H Tossavainen
Notes :
IIRC, They're also called walsh-hadamard transforms.
Basically like Fourier, but the basis functions are squarewaves with different sequencies.
I did this for a transform data compression study a while back.
Here's some code to do a walsh hadamard transform on long ints in-place (you need to divide by n to get transform) the order is bit-reversed at output, IIRC.
The inverse transform is the same as the forward transform (expects bit-reversed input). i.e. x = 1/n * FWHT(FWHT(x)) (x is a vector)
Code :
void inline wht_bfly (long& a, long& b)
{
long tmp = a;
a += b;
b = tmp - b;
}
// just a integer log2
int inline l2 (long x)
{
int l2;
for (l2 = 0; x > 0; x >>=1)
{
++ l2;
}
return (l2);
}
////////////////////////////////////////////
// Fast in-place Walsh-Hadamard Transform //
////////////////////////////////////////////
void FWHT (std::vector& data)
{
const int log2 = l2 (data.size()) - 1;
for (int i = 0; i < log2; ++i)
{
for (int j = 0; j < (1 << log2); j += 1 << (i+1))
{
for (int k = 0; k < (1<<i); ++k)
{
wht_bfly (data [j + k], data [j + k + (1<<i)]);
}
}
}
}
17 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
FFT
References : Toth Laszlo
Linked file : rvfft.ps
Linked file : rvfft.cpp
Notes :
A paper (postscript) and some C++ source for 4 different fft algorithms, compiled by Toth Laszlo from the Hungarian Academy of Sciences Research Group on Artificial Intelligence.
Toth says: "I've found that Sorensen's split-radix algorithm was the fastest, so I use this since then (this means that you may as well delete the other routines in my source - if you believe my results)."
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Frequency response from biquad coefficients
Type : biquad
References : Posted by peter[AT]sonicreef[DOT]com
Notes :
Here is a formula for plotting the frequency response of a biquad filter. Depending on the coefficients that you have, you might have to use negative values for the b- coefficients.
Code :
//w = frequency (0 < w < PI)
//square(x) = x*x
y = 20*log((sqrt(square(a0*square(cos(w))-a0*square(sin(w))+a1*cos(w)+a2)+square(2*a0*cos(w)*sin(w)+a1*(sin(w))))/
sqrt(square(square(cos(w))- square(sin(w))+b1*cos(w)+b2)+square(2* cos(w)*sin(w)+b1*(sin(w))))));
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Java FFT
Type : FFT Analysis
References : Posted by Loreno Heer
Notes :
May not work correctly ;-)
Code :
// WTest.java
/*
Copyright (C) 2003 Loreno Heer, (helohe at bluewin dot ch)
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
public class WTest{
private static double[] sin(double step, int size){
double f = 0;
double[] ret = new double[size];
for(int i = 0; i < size; i++){
ret[i] = Math.sin(f);
f += step;
}
return ret;
}
private static double[] add(double[] a, double[] b){
double[] c = new double[a.length];
for(int i = 0; i < a.length; i++){
c[i] = a[i] + b[i];
}
return c;
}
private static double[] sub(double[] a, double[] b){
double[] c = new double[a.length];
for(int i = 0; i < a.length; i++){
c[i] = a[i] - b[i];
}
return c;
}
private static double[] add(double[] a, double b){
double[] c = new double[a.length];
for(int i = 0; i < a.length; i++){
c[i] = a[i] + b;
}
return c;
}
private static double[] cp(double[] a, int size){
double[] c = new double[size];
for(int i = 0; i < size; i++){
c[i] = a[i];
}
return c;
}
private static double[] mul(double[] a, double b){
double[] c = new double[a.length];
for(int i = 0; i < a.length; i++){
c[i] = a[i] * b;
}
return c;
}
private static void print(double[] value){
for(int i = 0; i < value.length; i++){
System.out.print(i + "," + value[i] + "\n");
}
System.out.println();
}
private static double abs(double[] a){
double c = 0;
for(int i = 0; i < a.length; i++){
c = ((c * i) + Math.abs(a[i])) / (i + 1);
}
return c;
}
private static double[] fft(double[] a, int min, int max, int step){
double[] ret = new double[(max - min) / step];
int i = 0;
for(int d = min; d < max; d = d + step){
double[] f = sin(fc(d), a.length);
double[] dif = sub(a, f);
ret[i] = 1 - abs(dif);
i++;
}
return ret;
}
private static double[] fft_log(double[] a){
double[] ret = new double[1551];
int i = 0;
for(double d = 0; d < 15.5; d = d + 0.01){
double[] f = sin(fc(Math.pow(2,d)), a.length);
double[] dif = sub(a, f);
ret[i] = Math.abs(1 - abs(dif));
i++;
}
return ret;
}
private static double fc(double d){
return d * Math.PI / res;
}
private static void print_log(double[] value){
for(int i = 0; i < value.length; i++){
System.out.print(Math.pow(2,((double)i/100d)) + "," + value[i] + "\n");
}
System.out.println();
}
public static void main(String[] args){
double[] f_0 = sin(fc(440), sample_length); // res / pi =>14005
//double[] f_1 = sin(.02, sample_length);
double[] f_2 = sin(fc(520), sample_length);
//double[] f_3 = sin(.25, sample_length);
//double[] f = add( add( add(f_0, f_1), f_2), f_3);
double[] f = add(f_0, f_2);
//print(f);
double[] d = cp(f,1000);
print_log(fft_log(d));
}
static double length = .2; // sec
static int res = 44000; // resoultion (pro sec)
static int sample_length = res; // resoultion
}
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Look ahead limiting
References : Posted by Wilfried Welti
Notes :
use add_value with all values which enter the look-ahead area,
and remove_value with all value which leave this area. to get
the maximum value in the look-ahead area, use get_max_value.
in the very beginning initialize the table with zeroes.
If you always want to know the maximum amplitude in
your look-ahead area, the thing becomes a sorting
problem. very primitive approach using a look-up table
Code :
void lookup_add(unsigned section, unsigned size, unsigned value)
{
if (section==value)
lookup[section]++;
else
{
size >>= 1;
if (value>section)
{
lookup[section]++;
lookup_add(section+size,size,value);
}
else
lookup_add(section-size,size,value);
}
}
void lookup_remove(unsigned section, unsigned size, unsigned value)
{
if (section==value)
lookup[section]--;
else
{
size >>= 1;
if (value>section)
{
lookup[section]--;
lookup_remove(section+size,size,value);
}
else
lookup_remove(section-size,size,value);
}
}
unsigned lookup_getmax(unsigned section, unsigned size)
{
unsigned max = lookup[section] ? section : 0;
size >>= 1;
if (size)
if (max)
{
max = lookup_getmax((section+size),size);
if (!max) max=section;
}
else
max = lookup_getmax((section-size),size);
return max;
}
void add_value(unsigned value)
{
lookup_add(LOOKUP_VALUES>>1, LOOKUP_VALUES>>1, value);
}
void remove_value(unsigned value)
{
lookup_remove(LOOKUP_VALUES>>1, LOOKUP_VALUES>>1, value);
}
unsigned get_max_value()
{
return lookup_getmax(LOOKUP_VALUES>>1, LOOKUP_VALUES>>1);
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
LPC analysis (autocorrelation + Levinson-Durbin recursion)
References : Posted by mail[AT]mutagene[DOT]net
Notes :
The autocorrelation function implements a warped autocorrelation, so that frequency resolution can be specified by the variable 'lambda'. Levinson-Durbin recursion calculates autoregression coefficients a and reflection coefficients (for lattice filter implementation) K. Comments for Levinson-Durbin function implement matlab version of the same function.
No optimizations.
Code :
//find the order-P autocorrelation array, R, for the sequence x of length L and warping of lambda
//wAutocorrelate(&pfSrc[stIndex],siglen,R,P,0);
wAutocorrelate(float * x, unsigned int L, float * R, unsigned int P, float lambda)
{
double * dl = new double [L];
double * Rt = new double [L];
double r1,r2,r1t;
R[0]=0;
Rt[0]=0;
r1=0;
r2=0;
r1t=0;
for(unsigned int k=0; k<L;k++)
{
Rt[0]+=double(x[k])*double(x[k]);
dl[k]=r1-double(lambda)*double(x[k]-r2);
r1 = x[k];
r2 = dl[k];
}
for(unsigned int i=1; i<=P; i++)
{
Rt[i]=0;
r1=0;
r2=0;
for(unsigned int k=0; k<L;k++)
{
Rt[i]+=double(dl[k])*double(x[k]);
r1t = dl[k];
dl[k]=r1-double(lambda)*double(r1t-r2);
r1 = r1t;
r2 = dl[k];
}
}
for(i=0; i<=P; i++)
R[i]=float(Rt[i]);
delete[] dl;
delete[] Rt;
}
// Calculate the Levinson-Durbin recursion for the autocorrelation sequence R of length P+1 and return the autocorrelation coefficients a and reflection coefficients K
LevinsonRecursion(unsigned int P, float *R, float *A, float *K)
{
double Am1[62];
if(R[0]==0.0) {
for(unsigned int i=1; i<=P; i++)
{
K[i]=0.0;
A[i]=0.0;
}}
else {
double km,Em1,Em;
unsigned int k,s,m;
for (k=0;k<=P;k++){
A[0]=0;
Am1[0]=0; }
A[0]=1;
Am1[0]=1;
km=0;
Em1=R[0];
for (m=1;m<=P;m++) //m=2:N+1
{
double err=0.0f; //err = 0;
for (k=1;k<=m-1;k++) //for k=2:m-1
err += Am1[k]*R[m-k]; // err = err + am1(k)*R(m-k+1);
km = (R[m]-err)/Em1; //km=(R(m)-err)/Em1;
K[m-1] = -float(km);
A[m]=(float)km; //am(m)=km;
for (k=1;k<=m-1;k++) //for k=2:m-1
A[k]=float(Am1[k]-km*Am1[m-k]); // am(k)=am1(k)-km*am1(m-k+1);
Em=(1-km*km)*Em1; //Em=(1-km*km)*Em1;
for(s=0;s<=P;s++) //for s=1:N+1
Am1[s] = A[s]; // am1(s) = am(s)
Em1 = Em; //Em1 = Em;
}
}
return 0;
}
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Magnitude and phase plot of arbitrary IIR function, up to 5th order
Type : magnitude and phase at any frequency
References : Posted by George Yohng
Notes :
Amplitude and phase calculation of IIR equation
run at sample rate "sampleRate" at frequency "F".
AMPLITUDE
-----------
cf_mag(F,sampleRate,
a0,a1,a2,a3,a4,a5,
b0,b1,b2,b3,b4,b5)
-----------
PHASE
-----------
cf_phi(F,sampleRate,
a0,a1,a2,a3,a4,a5,
b0,b1,b2,b3,b4,b5)
-----------
If you need a frequency diagram, draw a plot for
F=0...sampleRate/2
If you need amplitude in dB, use cf_lin2db(cf_mag(.......))
Set b0=-1 if you have such function:
y[n] = a0*x[n] + a1*x[n-1] + a2*x[n-2] + a3*x[n-3] + a4*x[n-4] + a5*x[n-5] +
+ b1*y[n-1] + b2*y[n-2] + b3*y[n-3] + b4*y[n-4] + b5*y[n-5];
Set b0=1 if you have such function:
y[n] = a0*x[n] + a1*x[n-1] + a2*x[n-2] + a3*x[n-3] + a4*x[n-4] + a5*x[n-5] +
- b1*y[n-1] - b2*y[n-2] - b3*y[n-3] - b4*y[n-4] - b5*y[n-5];
Do not try to reverse engineer these formulae - they don't give any sense
other than they are derived from transfer function, and they work. :)
Code :
/*
C file can be downloaded from
http://www.yohng.com/dsp/cfsmp.c
*/
#define C_PI 3.14159265358979323846264
double cf_mag(double f,double rate,
double a0,double a1,double a2,double a3,double a4,double a5,
double b0,double b1,double b2,double b3,double b4,double b5)
{
return
sqrt((a0*a0 + a1*a1 + a2*a2 + a3*a3 + a4*a4 + a5*a5 +
2*(a0*a1 + a1*a2 + a2*a3 + a3*a4 + a4*a5)*cos((2*f*C_PI)/rate) +
2*(a0*a2 + a1*a3 + a2*a4 + a3*a5)*cos((4*f*C_PI)/rate) +
2*a0*a3*cos((6*f*C_PI)/rate) + 2*a1*a4*cos((6*f*C_PI)/rate) +
2*a2*a5*cos((6*f*C_PI)/rate) + 2*a0*a4*cos((8*f*C_PI)/rate) +
2*a1*a5*cos((8*f*C_PI)/rate) + 2*a0*a5*cos((10*f*C_PI)/rate))/
(b0*b0 + b1*b1 + b2*b2 + b3*b3 + b4*b4 + b5*b5 +
2*(b0*b1 + b1*b2 + b2*b3 + b3*b4 + b4*b5)*cos((2*f*C_PI)/rate) +
2*(b0*b2 + b1*b3 + b2*b4 + b3*b5)*cos((4*f*C_PI)/rate) +
2*b0*b3*cos((6*f*C_PI)/rate) + 2*b1*b4*cos((6*f*C_PI)/rate) +
2*b2*b5*cos((6*f*C_PI)/rate) + 2*b0*b4*cos((8*f*C_PI)/rate) +
2*b1*b5*cos((8*f*C_PI)/rate) + 2*b0*b5*cos((10*f*C_PI)/rate)));
}
double cf_phi(double f,double rate,
double a0,double a1,double a2,double a3,double a4,double a5,
double b0,double b1,double b2,double b3,double b4,double b5)
{
atan2((a0*b0 + a1*b1 + a2*b2 + a3*b3 + a4*b4 + a5*b5 +
(a0*b1 + a1*(b0 + b2) + a2*(b1 + b3) + a5*b4 + a3*(b2 + b4) +
a4*(b3 + b5))*cos((2*f*C_PI)/rate) +
((a0 + a4)*b2 + (a1 + a5)*b3 + a2*(b0 + b4) +
a3*(b1 + b5))*cos((4*f*C_PI)/rate) + a3*b0*cos((6*f*C_PI)/rate) +
a4*b1*cos((6*f*C_PI)/rate) + a5*b2*cos((6*f*C_PI)/rate) +
a0*b3*cos((6*f*C_PI)/rate) + a1*b4*cos((6*f*C_PI)/rate) +
a2*b5*cos((6*f*C_PI)/rate) + a4*b0*cos((8*f*C_PI)/rate) +
a5*b1*cos((8*f*C_PI)/rate) + a0*b4*cos((8*f*C_PI)/rate) +
a1*b5*cos((8*f*C_PI)/rate) +
(a5*b0 + a0*b5)*cos((10*f*C_PI)/rate))/
(b0*b0 + b1*b1 + b2*b2 + b3*b3 + b4*b4 + b5*b5 +
2*((b0*b1 + b1*b2 + b3*(b2 + b4) + b4*b5)*cos((2*f*C_PI)/rate) +
(b2*(b0 + b4) + b3*(b1 + b5))*cos((4*f*C_PI)/rate) +
(b0*b3 + b1*b4 + b2*b5)*cos((6*f*C_PI)/rate) +
(b0*b4 + b1*b5)*cos((8*f*C_PI)/rate) +
b0*b5*cos((10*f*C_PI)/rate))),
((a1*b0 + a3*b0 + a5*b0 - a0*b1 + a2*b1 + a4*b1 - a1*b2 +
a3*b2 + a5*b2 - a0*b3 - a2*b3 + a4*b3 -
a1*b4 - a3*b4 + a5*b4 - a0*b5 - a2*b5 - a4*b5 +
2*(a3*b1 + a5*b1 - a0*b2 + a4*(b0 + b2) - a1*b3 + a5*b3 +
a2*(b0 - b4) - a0*b4 - a1*b5 - a3*b5)*cos((2*f*C_PI)/rate) +
2*(a3*b0 + a4*b1 + a5*(b0 + b2) - a0*b3 - a1*b4 - a0*b5 - a2*b5)*
cos((4*f*C_PI)/rate) + 2*a4*b0*cos((6*f*C_PI)/rate) +
2*a5*b1*cos((6*f*C_PI)/rate) - 2*a0*b4*cos((6*f*C_PI)/rate) -
2*a1*b5*cos((6*f*C_PI)/rate) + 2*a5*b0*cos((8*f*C_PI)/rate) -
2*a0*b5*cos((8*f*C_PI)/rate))*sin((2*f*C_PI)/rate))/
(b0*b0 + b1*b1 + b2*b2 + b3*b3 + b4*b4 + b5*b5 +
2*(b0*b1 + b1*b2 + b2*b3 + b3*b4 + b4*b5)*cos((2*f*C_PI)/rate) +
2*(b0*b2 + b1*b3 + b2*b4 + b3*b5)*cos((4*f*C_PI)/rate) +
2*b0*b3*cos((6*f*C_PI)/rate) + 2*b1*b4*cos((6*f*C_PI)/rate) +
2*b2*b5*cos((6*f*C_PI)/rate) + 2*b0*b4*cos((8*f*C_PI)/rate) +
2*b1*b5*cos((8*f*C_PI)/rate) + 2*b0*b5*cos((10*f*C_PI)/rate)));
}
double cf_lin2db(double lin)
{
if (lin<9e-51) return -1000; /* prevent invalid operation */
return 20*log10(lin);
}
9 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Measuring interpollation noise
References : Posted by Jon Watte
Notes :
You can easily estimate the error by evaluating the actual function and
evaluating your interpolator at each of the mid-points between your
samples. The absolute difference between these values, over the absolute
value of the "correct" value, is your relative error. log10 of your relative
error times 20 is an estimate of your quantization noise in dB. Example:
You have a table for every 0.5 "index units". The value at index unit 72.0
is 0.995 and the value at index unit 72.5 is 0.999. The interpolated value
at index 72.25 is 0.997. Suppose the actual function value at that point was
0.998; you would have an error of 0.001 which is a relative error of 0.001002004..
log10(error) is about -2.99913, which times 20 is about -59.98. Thus, that's
your quantization noise at that position in the table. Repeat for each pair of
samples in the table.
Note: I said "quantization noise" not "aliasing noise". The aliasing noise will,
as far as I know, only happen when you start up-sampling without band-limiting
and get frequency aliasing (wrap-around), and thus is mostly independent of
what specific interpolation mechanism you're using.
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
QFT and DQFT (double precision) classes
References : Posted by Joshua Scholar
Linked file : qft.tar_1.gz
Notes :
Since it's a Visual C++ project (though it has relatively portable C++) I
guess the main audience are PC users. As such I'm including a zip file.
Some PC users wouldn't know what to do with a tgz file.
The QFT and DQFT (double precision) classes supply the following functions:
1. Real valued FFT and inverse FFT functions. Note that separate arraysare used
for real and imaginary component of the resulting spectrum.
2. Decomposition of a spectrum into a separate spectrum of the evensamples
and a spectrum of the odd samples. This can be useful for buildingfilter banks.
3. Reconstituting a spectrum from separate spectrums of the even samples
and odd samples. This can be useful for building filter banks.
4. A discrete Sin transform (a QFT decomposes an FFT into a DST and DCT).
5. A discrete Cos transfrom.
6. Since a QFT does it's last stage calculating from the outside in thelast part
can be left unpacked and only calculated as needed in the case wherethe entire
spectrum isn't needed (I used this for calculating correlations andconvolutions
where I only needed half of the results).
ReverseNoUnpack()
UnpackStep()
and NegUnpackStep()
implement this functionality
NOTE Reverse() normalizes its results (divides by one half the blocklength), but
ReverseNoUnpack() does not.
7. Also if you only want the first half of the results you can call ReverseHalf()
NOTE Reverse() normalizes its results (divides by one half the blocklength), but
ReverseHalf() does not.
8. QFT is less numerically stable than regular FFTs. With singleprecision calculations,
a block length of 2^15 brings the accuracy down to being barelyaccurate enough.
At that size, single precision calculations tested sound files wouldoccasionally have
a sample off by 2, and a couple off by 1 per block. Full volume whitenoise would generate
a few samples off by as much as 6 per block at the end, beginning and middle.
No matter what the inputs the errors are always at the same positions in the block.
There some sort of cancelation that gets more delicate as the block size gets bigger.
For the sake of doing convolutions and the like where the forward transform is
done only once for one of the inputs, I created a AccurateForward() function.
It uses a regular FFT algorithm for blocks larger than 2^12, and decomposes into even and
odd FFTs recursively.
In any case you can always use the double precision routines to get more accuracy.
DQFT even has routines that take floats as inputs and return double precision
spectrum outputs.
As for portability:
1. The files qft.cpp and dqft.cpp start with defines:
#define _USE_ASM
If you comment those define out, then what's left is C++ with no assembly language.
2. There is unnecessary windows specific code in "criticalSection.h"
I used a critical section because objects are not reentrant (each object has
permanent scratch pad memory), but obviously critical sections are operating
system specific. In any case that code can easily be taken out.
If you look at my code and see that there's an a test built in the examples
that makes sure that the results are in the ballpark of being right. It
wasn't that I expected the answers to be far off, it was that I uncommenting
the "no assembly language" versions of some routines and I wanted to make
sure that they weren't broken.
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Simple peak follower
Type : amplitude analysis
References : Posted by Phil Burk
Notes :
This simple peak follower will give track the peaks of a signal. It will rise rapidly when the input is rising, and then decay exponentially when the input drops. It can be used to drive VU meters, or used in an automatic gain control circuit.
Code :
// halfLife = time in seconds for output to decay to half value after an impulse
static float output = 0.0;
float scalar = pow( 0.5, 1.0/(halfLife * sampleRate)));
if( input < 0.0 )
input = -input; /* Absolute value. */
if ( input >= output )
{
/* When we hit a peak, ride the peak to the top. */
output = input;
}
else
{
/* Exponential decay of output when signal is low. */
output = output * scalar;
/*
** When current gets close to 0.0, set current to 0.0 to prevent FP underflow
** which can cause a severe performance degradation due to a flood
** of interrupts.
*/
if( output < VERY_SMALL_FLOAT ) output = 0.0;
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
tone detection with Goertzel
Type : Goertzel
References : Posted by espenr[AT]ii[DOT]uib[DOT]no
Linked file : http://www.ii.uib.no/~espenr/tonedetect.zip
Notes :
Goertzel is basically DFT of parts of a spectrum not the total spectrum as you normally do with FFT. So if you just want to check out the power for some frequencies this could be better. Is good for DTFM detection I've heard.
The WNk isn't calculated 100% correctly, but it seems to work so ;) Yeah and the code is C++ so you might have to do some small adjustment to compile it as C.
Code :
/** Tone detect by Goertzel algorithm
*
* This program basically searches for tones (sines) in a sample and reports the different dB it finds for
* different frequencies. Can easily be extended with some thresholding to report true/false on detection.
* I'm far from certain goertzel it implemented 100% correct, but it works :)
*
* Hint, the SAMPLERATE, BUFFERSIZE, FREQUENCY, NOISE and SIGNALVOLUME all affects the outcome of the reported dB. Tweak
* em to find the settings best for your application. Also, seems to be pretty sensitive to noise (whitenoise anyway) which
* is a bit sad. Also I don't know if the goertzel really likes float values for the frequency ... And using 44100 as
* samplerate for detecting 6000 Hz tone is kinda silly I know :)
*
* Written by: Espen Riskedal, espenr@ii.uib.no, july-2002
*/
#include <iostream>
#include <cmath>
#include <cstdlib>
using std::rand;
// math stuff
using std::cos;
using std::abs;
using std::exp;
using std::log10;
// iostream stuff
using std::cout;
using std::endl;
#define PI 3.14159265358979323844
// change the defines if you want to
#define SAMPLERATE 44100
#define BUFFERSIZE 8820
#define FREQUENCY 6000
#define NOISE 0.05
#define SIGNALVOLUME 0.8
/** The Goertzel algorithm computes the k-th DFT coefficient of the input signal using a second-order filter.
* http://ptolemy.eecs.berkeley.edu/papers/96/dtmf_ict/www/node3.html.
* Basiclly it just does a DFT of the frequency we want to check, and none of the others (FFT calculates for all frequencies).
*/
float goertzel(float *x, int N, float frequency, int samplerate) {
float Skn, Skn1, Skn2;
Skn = Skn1 = Skn2 = 0;
for (int i=0; i<N; i++) {
Skn2 = Skn1;
Skn1 = Skn;
Skn = 2*cos(2*PI*frequency/samplerate)*Skn1 - Skn2 + x[i];
}
float WNk = exp(-2*PI*frequency/samplerate); // this one ignores complex stuff
//float WNk = exp(-2*j*PI*k/N);
return (Skn - WNk*Skn1);
}
/** Generates a tone of the specified frequency
* Gotten from: http://groups.google.com/groups?hl=en&lr=&ie=UTF-8&oe=UTF-8&safe=off&selm=3c641e%243jn%40uicsl.csl.uiuc.edu
*/
float *makeTone(int samplerate, float frequency, int length, float gain=1.0) {
//y(n) = 2 * cos(A) * y(n-1) - y(n-2)
//A= (frequency of interest) * 2 * PI / (sampling frequency)
//A is in radians.
// frequency of interest MUST be <= 1/2 the sampling frequency.
float *tone = new float[length];
float A = frequency*2*PI/samplerate;
for (int i=0; i<length; i++) {
if (i > 1) tone[i]= 2*cos(A)*tone[i-1] - tone[i-2];
else if (i > 0) tone[i] = 2*cos(A)*tone[i-1] - (cos(A));
else tone[i] = 2*cos(A)*cos(A) - cos(2*A);
}
for (int i=0; i<length; i++) tone[i] = tone[i]*gain;
return tone;
}
/** adds whitenoise to a sample */
void *addNoise(float *sample, int length, float gain=1.0) {
for (int i=0; i<length; i++) sample[i] += (2*(rand()/(float)RAND_MAX)-1)*gain;
}
/** returns the signal power/dB */
float power(float value) {
return 20*log10(abs(value));
}
int main(int argc, const char* argv) {
cout << "Samplerate: " << SAMPLERATE << "Hz\n";
cout << "Buffersize: " << BUFFERSIZE << " samples\n";
cout << "Correct frequency is: " << FREQUENCY << "Hz\n";
cout << " - signal volume: " << SIGNALVOLUME*100 << "%\n";
cout << " - white noise: " << NOISE*100 << "%\n";
float *tone = makeTone(SAMPLERATE, FREQUENCY, BUFFERSIZE, SIGNALVOLUME);
addNoise(tone, BUFFERSIZE,NOISE);
int stepsize = FREQUENCY/5;
for (int i=0; i<10; i++) {
int freq = stepsize*i;
cout << "Trying freq: " << freq << "Hz -> dB: " << power(goertzel(tone, BUFFERSIZE, freq, SAMPLERATE)) << endl;
}
delete tone;
return 0;
}
10 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Tone detection with Goertzel (x86 ASM)
Type : Tone detection with Goertzel in x86 assembly
References : Posted by Christian[AT]savioursofsoul[DOT]de
Notes :
This is an "assemblified" version of the Goertzel Tone Detector. It is about 2 times faster than the original code.
The code has been tested and it works fine.
Hope you can use it. I'm gonna try to build a Tuner (as VST-Plugin). I hope, that this will work :-\ If anyone is intrested, please let me know.
Christian
Code :
function Goertzel_x87(Buffer :Psingle; BLength:Integer; frequency: Single; samplerate: Single):Single;
asm
mov ecx,BLength
mov eax,Buffer
fld x2
fldpi
fmulp
fmul frequency
fdiv samplerate
fld st(0)
fcos
fld x2
fmulp
fxch st(1)
fldz
fsub st(0),st(1)
fstp st(1)
fldl2e
fmul
fld st(0)
frndint
fsub st(1),st(0)
fxch st(1)
f2xm1
fld1
fadd
fscale
fstp st(1)
fldz
fldz
fldz
@loopStart:
fxch st(1)
fxch st(2)
fstp st(0)
fld st(3)
fmul st(0),st(1)
fsub st(0),st(2)
fld [eax].Single
faddp
add eax,4
loop @loopStart
@loopEnd:
fxch st(3)
fmulp st(2), st(0)
fsub st(0),st(1)
fstp result
ffree st(2)
ffree st(1)
ffree st(0)
end;
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Vintage VU meters tutorial
References : Posted by neolit123 at gmail dot com
Notes :
Here is a short tutorial about vintage-styled VU meters:
http://neolit123.blogspot.com/2009/03/designing-analog-vu-meter-in-dsp.html
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
