TOPIC: Re: Digital crossover for my WMTMW speakers

clf(1)

//Definitions
ntaps = 500
fs = 44100
freqvect = [0:ntaps-1]*fs/2/ntaps

//Import measured data and the desired response (target)
//Measure data is assumed to be three columns of data: frequency (Hz), magnitude (dB), phase (rad).

//For the target however, I choose to have a linear scale for the magnitude, since I'd like zero for the stopband (which gets -inf in dB)

measureddata = fscanfMat('/path/measured.txt');
targetdata = fscanfMat('/path/target.txt');

//For calculation of the "ideal" filter - the imported target data is interpolated to the same frequency vector as the measured data

f=measureddata(:,1) //The frequency vector

//Reform the measured data into complex rectangular form
magn = measureddata(:,2)
phase = measureddata(:,3)
Repart = 10^(magn/20).*cos(phase*2*%pi/360)
Impart = 10^(magn/20).*sin(phase*2*%pi/360)
measured_complex = [Repart+Impart*%i];

//Interpolate the target function for the same frequency
df = splin(targetdata(:,1),targetdata(:,2))
dp = splin(targetdata(:,1),targetdata(:,3))
magn = interp(f,targetdata(:,1),targetdata(:,2),df)
phase = interp(f, targetdata(:,1),targetdata(:,3),dp)
Repart = magn.*cos(phase*2*%pi/360)
Impart = magn.*sin(phase*2*%pi/360)
target_complex = [Repart+Impart*%i];

//Now, create the "ideal" filter by dividing the target with the measured data
ideal_complex = target_complex./measured_complex

//The ideal filter is truncated in the time domain by resampling in the frequency domain. The interp function require real numbers, so the complex ideal filter must first be separated into real and imaginary parts
Re=real(ideal_complex)
Im=imag(ideal_complex)

dr = splin(f,Re)
di = splin(f,Im)
Re_trunc = interp(freqvect,f,Re,dr)
Im_trunc = interp(freqvect,f,Im,di)

//Combine the resampled real vectors into the new complex filter response vector
filter_trunc = Re_trunc+Im_trunc*%i

//Do ifft on the filter response and shift the result
filter_trunc_ifft = (fftshift(ifft(filter_trunc)))

//Apply the window to the truncated filter (in this case a Hamming window)
wind_filter=filter_trunc_ifft.*window('kr',ntaps,5)

//Transform back to get the frequency and phase response of the filter
filter_response = fft(wind_filter)

//Calculating the magnitude response of the filter
filter_magn = abs(filter_response(1:ntaps))

//Calculating the magnitude response of the "ideal" filter
ideal_response = abs(ideal_complex)

// If there are any zero elements the logarithm can't be calculated, so zeros are replaced by a very small number
x=find (ideal_response<=0.) ; ideal_response(x)=0.00000001
x=find (targetdata(:,2)<=0.) ; targetdata(x,2)=0.00000001

//Calculate the resulting response after filtering by convolution in the time domain
Result = fft(convol((wind_filter),fftshift(ifft(measured_complex))))

//Plot the result
scf(1)
subplot(2,3,1)
plot2d ('ln',measureddata(:,1),measureddata(:,2),rect=[20.,-100.,20000.,10.]);title("Imported measured response",'fontsize',3)
subplot(2,3,2)
plot2d ('ln',targetdata(:,1),20*log10(targetdata(:,2)),rect=[20.,-100.,20000.,10.]);title("Imported target response",'fontsize',3)
subplot(2,3,3)
plot2d ('ln',f,20*log10(ideal_response),rect=[20.,-100.,20000.,10.]);title("Ideal filter response",'fontsize',3)
subplot(2,3,4)
plot2d(wind_filter);title("Truncated filter",'fontsize',3)
subplot(2,3,5)
plot2d ('ln',freqvect,20*log10(filter_magn),rect=[20.,-100.,20000.,10.]);title("Resulting filter response",'fontsize',3)
subplot(2,3,6)
plot2d ('ln',f,20*log10(abs(Result(1:length(f)))),rect=[20.,-100.,20000.,10.]);title("Result of convolution measured data with filter",'fontsize',3)

f2 = [0:length(f)+ntaps-2]*fs/2/(length(f)+ntaps-1)

f

targetdata = [measureddata(:,1),(measureddata(:,1)<=50)*1.+(measureddata(:,1)>50)*0.,measureddata(:,1)*0.]

//Definitions for Keele-Horbach filter
WLSpace = 0.4 //Assumed wavelength spacing of speakers
SpdSnd = 343.2 //Speed of sound (m/s)
WfDist = 0.5 //Distance (m) between centers of outer spaced speakers (woofers)
MidDist = 0.15 //Distance (m) between centers of inner spaced speakers (midranges)
StepRatio = WfDist/MidDist
LPCritfreq = WLSpace*SpdSnd/WfDist
BPCritfreq = StepRatio*LPCritfreq

LPresponse = freqvect*0+(freqvect<=LPCritfreq)*1+((freqvect>LPCritfreq)&(freqvect<=LPCritfreq*StepRatio))*0.5.*abs((2*cos(1/3*%pi*freqvect/LPCritfreq/StepRatio)-1)./(cos(1/3*%pi*freqvect/LPCritfreq)-cos(1/3*%pi*freqvect/LPCritfreq./StepRatio)+1*(freqvect==0)))+(freqvect>LPCritfreq*StepRatio)*0

BPresponse = freqvect*0+(freqvect<=LPCritfreq*StepRatio).*(1-LPresponse)+((freqvect>BPCritfreq)&(freqvect<=3*BPCritfreq)).*(1-cos(%pi/3*freqvect/BPCritfreq)+1*(freqvect==0))^(-1)/2//-((freqvect>3*BPCritfreq)&(freqvect<=6*BPCritfreq)).*(0.25*freqvect/BPCritfreq/3-0.5) //This last term is for smooting the transition between midrange and tweeter, see fig. 13 in the 2nd paper by Horbach & Keele

HPresponse = 1-LPresponse-BPresponse

//Truncate in the time domain
filter_trunc_ifft=fftshift(ifft(ideal_complex))
filter_trunc_ifft((length(f)+ntaps)/2:length(f)-1) = []
filter_trunc_ifft(1:(length(f)-ntaps)/2)=[]

//Apply the window to the truncated filter (in this case a Hamming window)
wind_filter=filter_trunc_ifft.*window('kr',ntaps,5)'

atan(imag(Result),real(Result)

Phase=atan(imag(Result)./real(Result))