DSPRelated.com
Code
   Code Snippets    Submitted by hordia

Block processing adaptative notch filter

Hernán Ordiales ● December 17, 2010 ● Coded in C++ 
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

Hernán Ordiales ● November 19, 2010 ● Coded in Python 
#!/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

Hernán Ordiales ● November 9, 2010 ● Coded in Matlab 
%
%   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