Welcome, Guest
Username: Password: Remember me

TOPIC: Re: Digital crossover for my WMTMW speakers

Re: Digital crossover for my WMTMW speakers 7 years 8 months ago #4617

  • Crumboo
  • Crumboo's Avatar
  • Offline
  • Expert Boarder
  • Posts: 144
  • Thank you received: 14
  • Karma: 13
Hi again,

I've been doing some more work in Scilab and have some preliminary results. In short, the script works like this:
  • Import measured data and target response (text files with 3 columns: frequency, magnitude and phase)
  • Create the "ideal" filter by dividing the target response by the measured data
  • Truncate the ideal filter by resampling in the frequency domain
  • Do IFFT and apply the window
  • FFT back in order to evaluate the filter response

Here is one example: I used measured near-field response for my midrange and the target response from Keele/Horbach bandpass.


Here are the results for filters using 500, 1000 and 2000 taps:




There is still work to do; for example I haven't yet looked at the phase response. But it's a start. If you want I can share the code.

Best regards. :)
Last Edit: 7 years 8 months ago by Crumboo.
The administrator has disabled public write access.

Re: Digital crossover for my WMTMW speakers 7 years 8 months ago #4618

  • Jakster
  • Jakster's Avatar
  • Offline
  • Junior Boarder
  • Posts: 24
  • Thank you received: 4
  • Karma: 3
Awesome!!
If you want I can share the code.

yes please.. I´ve spent way to much time in scilab trying to get functions to work. I wanted subfunctions I could call from a main program. Haven´t made it work yet.

Filters look good - have you tried verifying by convoluting speakerrespone with correction filter and comparing it with you desired response (all in time domain)??

Also as Curryman pointed out earlier, the longest filters are used when operating at low frequencies. Have you tried working with the bas response (with at cutoff at around 30Hz) - I´d be interested in seeing how big filters are necessary for a fast cut off.

Great progress - keep it up.
The administrator has disabled public write access.

Re: Digital crossover for my WMTMW speakers 7 years 8 months ago #4619

  • Crumboo
  • Crumboo's Avatar
  • Offline
  • Expert Boarder
  • Posts: 144
  • Thank you received: 14
  • Karma: 13
Filters look good - have you tried verifying by convoluting speakerrespone with correction filter and comparing it with you desired response (all in time domain)??

Yes, here it comes - looks good (same filter as previous post, 500 taps BP):


I'll have a look on a low frequency filter soon - or if you want to do it yourself; here is the code:
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)

B)
The administrator has disabled public write access.
The following user(s) said Thank You: Jakster

Re: Digital crossover for my WMTMW speakers 7 years 8 months ago #4620

  • Jakster
  • Jakster's Avatar
  • Offline
  • Junior Boarder
  • Posts: 24
  • Thank you received: 4
  • Karma: 3
somebody´s been busy! :woohoo: my parental leave was spent carrying around a screaming baby... I would rather have spent it like you do :blink:

That file looks really good - It´s a shame I can´t give you more karma (once every six hours). I like the fact that you already have incorporated interpolation to match your target and measurement files (if they are different lengths). I like your use of the kaiser window for reducing out of band sidelobes.

I went (briefly) through your code and can´t see any errors. Would it be too much to ask if you´d upload your target filter text file?
The administrator has disabled public write access.

Re: Digital crossover for my WMTMW speakers 7 years 8 months ago #4621

  • Crumboo
  • Crumboo's Avatar
  • Offline
  • Expert Boarder
  • Posts: 144
  • Thank you received: 14
  • Karma: 13
Hi again - and thank you for your response! This work is really fun, and I'm glad the little one spends so much time sleeping ;)

I have to do a correction for the last plot. Since the length of the result from the convoluted response is longer than the input, one should define a new frequency vector for that plot (the problem gets very obvious when using a number of taps that's in the same order of magnitude as the number of samples in the data). To solve this, define the following frequency vector:
f2 = [0:length(f)+ntaps-2]*fs/2/(length(f)+ntaps-1)

and use it for the last plot instead of
f
.

I did some tests with a brickwall lowpass filter at 50 Hz. You can try for yourself by using this instead of the imported target:
targetdata = [measureddata(:,1),(measureddata(:,1)<=50)*1.+(measureddata(:,1)>50)*0.,measureddata(:,1)*0.]

Tried with 5000, 10000 and 30000 taps:

5000 taps:

10000 taps:

30000 taps:


Also, I attached my .txt files.

Good luck! :)


By the way, if one want to calculate the Keele-Horbach functions this script can be used (see their paper for info about the parameters):
//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
Attachments:
Last Edit: 7 years 8 months ago by Crumboo. Reason: Correcting the images
The administrator has disabled public write access.

Re: Digital crossover for my WMTMW speakers 7 years 8 months ago #4622

  • Crumboo
  • Crumboo's Avatar
  • Offline
  • Expert Boarder
  • Posts: 144
  • Thank you received: 14
  • Karma: 13
Hmm...found another thing that I believe is an error. Look at the resulting response after convolution in this plot (down to the right):


The irregularities in the speaker response between 2000 and 4000 Hz should have been corrected for by the filter, but that's obviously not the case here. Maybe there is a problem of truncating the "ideal" filters response by resampling in the frequency domain - I believed it should be equivalent to do a truncation in the time domain. Hope we can sort this out!
Last Edit: 7 years 8 months ago by Crumboo.
The administrator has disabled public write access.

Re: Digital crossover for my WMTMW speakers 7 years 8 months ago #4625

  • curryman
  • curryman's Avatar
  • Offline
  • Platinum Boarder
  • Posts: 791
  • Thank you received: 181
  • Karma: 100
Great progress! I'll try to find some time this evening to follow your work. In my company we also have tools to develop FIR filters. I'll try to check your work for consistency with the results of our software tools.

Again: great work!
The administrator has disabled public write access.

Re: Digital crossover for my WMTMW speakers 7 years 8 months ago #4626

  • Crumboo
  • Crumboo's Avatar
  • Offline
  • Expert Boarder
  • Posts: 144
  • Thank you received: 14
  • Karma: 13
@curryman: That would be very useful, thanks!


I think I've fixed the probelm - I changed the truncation to be in the time domain. Now it looks better: :)




Use this code for the truncation and windowing instead:
//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)'

So, the previously plot I posted are wrong I'm afraid ;)
The administrator has disabled public write access.

Re: Digital crossover for my WMTMW speakers 7 years 8 months ago #4629

  • Crumboo
  • Crumboo's Avatar
  • Offline
  • Expert Boarder
  • Posts: 144
  • Thank you received: 14
  • Karma: 13
I just keep on spamming here... ;)

Did another try with a brick wall LP at 50 Hz, 6000 taps. Looks like this:




Another thing: the imported data is assumed to be equally spaced in frequency (linear) and cover the whole range 0 Hz to fs/2. I guess using data that is logarithmically spaced will give errors... You can use for example Speaker Workshop for export linearly spaced data.
The administrator has disabled public write access.

Re: Digital crossover for my WMTMW speakers 7 years 8 months ago #4630

  • curryman
  • curryman's Avatar
  • Offline
  • Platinum Boarder
  • Posts: 791
  • Thank you received: 181
  • Karma: 100
Wow, lot`s of work.

I got everything working in Scilab and could follow all your tests. I also imported the measured response and the target filter to our companies software (HEAD acoustics ArtemiS). Took me some time to get the format of the files right for the import to ArtemiS, but now the data looks consistent to your graphs:




Next time i'll calculate the resulting FIR Filter, but for today I give up (too tired ;) )

This message has an attachment image.
Please log in or register to see it.

The administrator has disabled public write access.

Re: Digital crossover for my WMTMW speakers 7 years 8 months ago #4635

  • curryman
  • curryman's Avatar
  • Offline
  • Platinum Boarder
  • Posts: 791
  • Thank you received: 181
  • Karma: 100
Hi Crumboo, hi all

Today I calculated the target FIR Filter with ArtemiS -> assuming that our companies software is well tested (hope so ;)), crumboos Scilab results look identical and thus should be correct :)

Here's the resulting filter for 512 taps:

This message has an attachment image.
Please log in or register to see it.

The administrator has disabled public write access.
The following user(s) said Thank You: Crumboo

Re: Digital crossover for my WMTMW speakers 7 years 8 months ago #4636

  • Crumboo
  • Crumboo's Avatar
  • Offline
  • Expert Boarder
  • Posts: 144
  • Thank you received: 14
  • Karma: 13
Thank you curryman, that was good news! :)
The administrator has disabled public write access.

Re: Digital crossover for my WMTMW speakers 7 years 8 months ago #4647

  • Crumboo
  • Crumboo's Avatar
  • Offline
  • Expert Boarder
  • Posts: 144
  • Thank you received: 14
  • Karma: 13
I have been thinking about phase - I get this weird looking phase response from the result of the convolution (the phase was calculated as
atan(imag(Result),real(Result)
):



Any clues on this? As the phase of the target is zero for all frequencies: should one expect the phase of the result also to be zero - or some delay? How is the delay of the filter calculated (is it maybe as easy as the number of taps divided by the sample frequency)? Is it possible to plot the filter phase in the software you used curryman?
The administrator has disabled public write access.

Re: Digital crossover for my WMTMW speakers 7 years 8 months ago #4650

  • curryman
  • curryman's Avatar
  • Offline
  • Platinum Boarder
  • Posts: 791
  • Thank you received: 181
  • Karma: 100
sure ;)

Your measured response incl. phase:



Your target response incl. phase:



The resulting filter incl. phase:



The phase of the resulting filter looks a bit weird below 300Hz and above 6kHz, though not that weird as your plot :unsure: . This is caused by the phase response of the target that jumps between -180 and 0 in this frequency range (after interpolation, original phase response of uninterpolated target seems OK -> we should check the Scilab code here).

In Matlab phase response can be shown with unwrap(angle(complex vector)), not sure if this is available in Scilab (can check this this evening).

This message has attachments images.
Please log in or register to see it.

The administrator has disabled public write access.
The following user(s) said Thank You: Crumboo

Re: Digital crossover for my WMTMW speakers 7 years 8 months ago #4652

  • Crumboo
  • Crumboo's Avatar
  • Offline
  • Expert Boarder
  • Posts: 144
  • Thank you received: 14
  • Karma: 13
Great curryman! :)

Hm, I tried using the following expression for calculating the phase instead
Phase=atan(imag(Result)./real(Result))

and it seems to give a better result.

Measured phase (-pi to pi):

Phase of the filter:

Phase of the result, which is the impulse response of the data convoluted with the filter:


So, the result seems to be quite close to zero in the pass band... The problem in the stop band I believe, is that the complex target has zero as both real and imaginary parts. The calculation of phase will then give errors because of the division by zero.

Is it possible to do the convolution in your software curryman and check for the resulting phase response?
The administrator has disabled public write access.
Moderators: devteam