Code Snippets Submitted by hordia

## Block processing adaptative notch filter

```
bool Do(const Audio& input, const Audio& reference, Audio& output)
{
const unsigned size = input.GetSize();
const DataArray& in = input.GetBuffer();
const DataArray& ref = reference.GetBuffer();
DataArray& out = output.GetBuffer();
const unsigned M = mConfig.GetFilterLength();
const TData mu = mConfig.GetStep();
const unsigned L = M-1;
std::vector<TData> &wn = _filterCoef;
TData r[M];
for (unsigned k=0;k<size;k++)
{
TData acum = 0.;
for (unsigned n=0, i=k;n<M;n++,i++)
{
if (i<L)
r[n] = _oldRef[i];
else
r[n] = ref[i-L];
acum += wn[n]*r[n]; // interference estimation
}
if (k<L)
out[k] = _oldIn[k]-acum; // ECG signal estimation
else
out[k] = in[k-L]-acum; // ECG signal estimation
for (unsigned n=0;n<M;n++) wn[n]+=mu*r[n]*out[k]; // new filter taps
}
for (unsigned k=0;k<L;k++)
{
_oldRef[k] = ref[size-L+k];
_oldIn[k] = in[size-L+k];
}
return true;
}
```

## Fast audio convolution

```
#!/usr/bin/python
# -*- coding: utf-8 -*-
#
# Copyright (C) 2006 by Hernán Ordiales
# <audiocode@uint8.com.ar>
#
# 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.
#
# Fast convolution of an audio signal with the impulse response of a system modeled as LTI (Linear Time Invariant) one.
# Use: fast_conv.py source.wav impulse_response.wav output.wav
from wav_array import *
from FFT import fft, inverse_fft
from pylab import size
import sys
def nextpow2( L ):
N = 2
while N < L: N = N * 2
return N
def fast_conv_vect( x, h ):
# searches for the amount of points required to perform the FFT
L = size(h) + size(x) - 1 # linear convolution length
N = nextpow2( L )
# Note: N>=L is needed because the IDFT of the multiplication is the circular convolution and to match it to the
# common one, N>=L is required (where L=N1+N2-1;N1=length(x);N2=length(h))
# FFT(X,N) is the N points FFT, with zero padding if X has less than N points and truncated if has more.
H = fft( h, N ) # Fourier transform of the impulse
X = fft( x, N ) # Fourier transform of the input signal
Y = H * X # spectral multiplication
y = inverse_fft( Y ) # time domain again
return y
source = sys.argv[1]
impulse = sys.argv[2]
conv_output = sys.argv[3]
clip_factor = 1.01 # clip factor default value
[ h1, Fs1, h_bits ] = wavread( impulse ) # impulse response
[ x1, Fs2, x_bits ] = wavread( source ) # file to process
if Fs1 == Fs2 : # if sample rates are the same
print "Processing..."
y1 = fast_conv_vect( x1, h1 ).real # takes the real part to avoid a too small complex part (around e-18)
# audio normalization: if "y = y/max(y)" -> "values out of [-1,+1] range are clipped"
y1 = y1/( max(y1)*clip_factor ) # to avoid clipping
wavwrite( y1, Fs1, conv_output )
print "Output file:", conv_output
print "Convolution success."
else:
print "Error: files has different sample rate."
```

## Fast convolution

```
%
% Copyright (C) 2004 by Hernán Ordiales
% <audiocode@uint8.com.ar>
%
% 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.
%
function fast_conv( file, impulse, output )
% Fast convolution of an audio signal with the impulse response of a system modeled as LTI (Linear Time Invariant) one.
% USE: fast_conv( 'source.wav', 'impulse_response.wav', 'output.wav' );
h_stereo = true;
x_stereo = true;
clip_factor = 1.01; % clip factor default value
% impulse response:
[ h, Fs1 ] = wavread( impulse ); % 2 columns matrix, one per channel
h1 = h(:,1); % col 1 -> left channel
h_channels = size(h); h_channels = h_channels(2);
if( h_channels == 2 )
h2 = h(:,2); % col 2 -> right channel
else
h2 = h(:,1); % col 1 -> left channel
h_stereo = false;
end
% file to process:
[ x, Fs2 ] = wavread( file );
x1 = x(:,1); % col 1 -> left channel
x_channels = size(x); x_channels = x_channels(2);
if( x_channels == 2 )
x2 = x(:,2); % col 2 -> right channel
else
x2 = x(:,1); % col 1 -> left channel
x_stereo = false;
end
% if sample rates are the same
if( Fs1 == Fs2 )
% searches for the amount of points required to perform the FFT
L = length(h1) + length(x1) - 1; % linear convolution length
disp('Processing...');
% Use nextpow if you have it to replace below code -> N = 2^nextpow(L);
N = 2;
while( N<L ) % first 2 pow greater than L
N = N*2;
end
% Note: N>=L is needed because the IDFT of the multiplication is the circular convolution and to match it to the
% common one, N>=L is required (where L=N1+N2-1;N1=length(x);N2=length(h))
% FFT(X,N) is the N points FFT, with zero padding if X has less than N points and truncated if has more.
H1 = fft( h1, N ); % Fourier transform of the impulse
X1 = fft( x1, N ); % Fourier transform of the input signal
F_y1 = H1 .* X1; % spectral multiplication
y1 = ifft( F_y1 ); % time domain again
% audio normalization: if "y = y/max(y)" -> "values out of [-1,+1] range are clipped"
y1 = y1/( max(y1)*clip_factor ); % to avoid clipping
% if one of the files is stereo, process the other channel
if( h_stereo || x_stereo )
H2 = fft( h2, N ); % Fourier transform of the impulse
X2 = fft( x2, N ); % Fourier transform of the input signal
F_y2 = H2 .* X2; % spectral multiplication
y2 = ifft( F_y2 ); % time domain again
y2 = y2/( max(y2)*clip_factor ); % to avoid clipping
y = [ y1 y2 ]; % stereo output (two channels)
else
y = y1; % mono output
end
wavwrite( y, Fs1, output );
disp( 'Convolution success.' );
else
disp('Error: files has different sample rate.');
end
```

## Block processing adaptative notch filter

```
bool Do(const Audio& input, const Audio& reference, Audio& output)
{
const unsigned size = input.GetSize();
const DataArray& in = input.GetBuffer();
const DataArray& ref = reference.GetBuffer();
DataArray& out = output.GetBuffer();
const unsigned M = mConfig.GetFilterLength();
const TData mu = mConfig.GetStep();
const unsigned L = M-1;
std::vector<TData> &wn = _filterCoef;
TData r[M];
for (unsigned k=0;k<size;k++)
{
TData acum = 0.;
for (unsigned n=0, i=k;n<M;n++,i++)
{
if (i<L)
r[n] = _oldRef[i];
else
r[n] = ref[i-L];
acum += wn[n]*r[n]; // interference estimation
}
if (k<L)
out[k] = _oldIn[k]-acum; // ECG signal estimation
else
out[k] = in[k-L]-acum; // ECG signal estimation
for (unsigned n=0;n<M;n++) wn[n]+=mu*r[n]*out[k]; // new filter taps
}
for (unsigned k=0;k<L;k++)
{
_oldRef[k] = ref[size-L+k];
_oldIn[k] = in[size-L+k];
}
return true;
}
```

## Fast audio convolution

```
#!/usr/bin/python
# -*- coding: utf-8 -*-
#
# Copyright (C) 2006 by Hernán Ordiales
# <audiocode@uint8.com.ar>
#
# 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.
#
# Fast convolution of an audio signal with the impulse response of a system modeled as LTI (Linear Time Invariant) one.
# Use: fast_conv.py source.wav impulse_response.wav output.wav
from wav_array import *
from FFT import fft, inverse_fft
from pylab import size
import sys
def nextpow2( L ):
N = 2
while N < L: N = N * 2
return N
def fast_conv_vect( x, h ):
# searches for the amount of points required to perform the FFT
L = size(h) + size(x) - 1 # linear convolution length
N = nextpow2( L )
# Note: N>=L is needed because the IDFT of the multiplication is the circular convolution and to match it to the
# common one, N>=L is required (where L=N1+N2-1;N1=length(x);N2=length(h))
# FFT(X,N) is the N points FFT, with zero padding if X has less than N points and truncated if has more.
H = fft( h, N ) # Fourier transform of the impulse
X = fft( x, N ) # Fourier transform of the input signal
Y = H * X # spectral multiplication
y = inverse_fft( Y ) # time domain again
return y
source = sys.argv[1]
impulse = sys.argv[2]
conv_output = sys.argv[3]
clip_factor = 1.01 # clip factor default value
[ h1, Fs1, h_bits ] = wavread( impulse ) # impulse response
[ x1, Fs2, x_bits ] = wavread( source ) # file to process
if Fs1 == Fs2 : # if sample rates are the same
print "Processing..."
y1 = fast_conv_vect( x1, h1 ).real # takes the real part to avoid a too small complex part (around e-18)
# audio normalization: if "y = y/max(y)" -> "values out of [-1,+1] range are clipped"
y1 = y1/( max(y1)*clip_factor ) # to avoid clipping
wavwrite( y1, Fs1, conv_output )
print "Output file:", conv_output
print "Convolution success."
else:
print "Error: files has different sample rate."
```

## Fast convolution

```
%
% Copyright (C) 2004 by Hernán Ordiales
% <audiocode@uint8.com.ar>
%
% 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.
%
function fast_conv( file, impulse, output )
% Fast convolution of an audio signal with the impulse response of a system modeled as LTI (Linear Time Invariant) one.
% USE: fast_conv( 'source.wav', 'impulse_response.wav', 'output.wav' );
h_stereo = true;
x_stereo = true;
clip_factor = 1.01; % clip factor default value
% impulse response:
[ h, Fs1 ] = wavread( impulse ); % 2 columns matrix, one per channel
h1 = h(:,1); % col 1 -> left channel
h_channels = size(h); h_channels = h_channels(2);
if( h_channels == 2 )
h2 = h(:,2); % col 2 -> right channel
else
h2 = h(:,1); % col 1 -> left channel
h_stereo = false;
end
% file to process:
[ x, Fs2 ] = wavread( file );
x1 = x(:,1); % col 1 -> left channel
x_channels = size(x); x_channels = x_channels(2);
if( x_channels == 2 )
x2 = x(:,2); % col 2 -> right channel
else
x2 = x(:,1); % col 1 -> left channel
x_stereo = false;
end
% if sample rates are the same
if( Fs1 == Fs2 )
% searches for the amount of points required to perform the FFT
L = length(h1) + length(x1) - 1; % linear convolution length
disp('Processing...');
% Use nextpow if you have it to replace below code -> N = 2^nextpow(L);
N = 2;
while( N<L ) % first 2 pow greater than L
N = N*2;
end
% Note: N>=L is needed because the IDFT of the multiplication is the circular convolution and to match it to the
% common one, N>=L is required (where L=N1+N2-1;N1=length(x);N2=length(h))
% FFT(X,N) is the N points FFT, with zero padding if X has less than N points and truncated if has more.
H1 = fft( h1, N ); % Fourier transform of the impulse
X1 = fft( x1, N ); % Fourier transform of the input signal
F_y1 = H1 .* X1; % spectral multiplication
y1 = ifft( F_y1 ); % time domain again
% audio normalization: if "y = y/max(y)" -> "values out of [-1,+1] range are clipped"
y1 = y1/( max(y1)*clip_factor ); % to avoid clipping
% if one of the files is stereo, process the other channel
if( h_stereo || x_stereo )
H2 = fft( h2, N ); % Fourier transform of the impulse
X2 = fft( x2, N ); % Fourier transform of the input signal
F_y2 = H2 .* X2; % spectral multiplication
y2 = ifft( F_y2 ); % time domain again
y2 = y2/( max(y2)*clip_factor ); % to avoid clipping
y = [ y1 y2 ]; % stereo output (two channels)
else
y = y1; % mono output
end
wavwrite( y, Fs1, output );
disp( 'Convolution success.' );
else
disp('Error: files has different sample rate.');
end
```