Help on making MULTIRESOLUTION WAVELET ANALYSIS available ?

mige0 · 2010-04-19 9:57 pm

function [wtf_mag_norm] = wtf_mag_norm(wtf_mag, d_bug)

%--------------------------------------------------------
% wtf_mag_norm
%--------------------------------------------------------
% This file is part of Diy_Audio_WTF package, intended to get started on Audio Wavelet Transformation
% on free Octave software
% Copyright (C) 2010 Michael Gerstgrasser
%
% "wtf_mag_norm" is a function to normalize magnitude of wavelet analysis
%
% Input variables :
% wtf_mag is magnitude of wavelet analysis
% d_bug is for providing additional information in case of debug mode
%
% Output variables:
% wtf_mag_norm is normalized magnitude of wavelet analysis
%
%
% see also load_wav_ir.m, load_ascii_ir.m, load_ascii1_ir.m, fft_ir_raw, waveletdemo.m

% Contributing authors to the Diy_Audio_WTF package in order of code publishing date:
%
% Michael Gerstgrasser - Contact information: user 'mige0' at www.diyaudio.com
% Elias Pekonen - Contact information: user 'Elias' at www.diyaudio.com
% J. Lewalle: Prof. Jacques Lewalle
%
% check out on www.diyaudio.com :
% "Help on making MULTIRESOLUTION WAVELET ANALYSIS available ?"
% http://www.diyaudio.com/forums/soft...ltiresolution-wavelet-analysis-available.html
% "WTF!? Wavelet TransForm for audio measurements - What-is? and How-to?"
% http://www.diyaudio.com/forums/mult...et-transform-audio-measurements-what-how.html
% "Horn Honk $$ WANTED $$"
% http://www.diyaudio.com/forums/multi-way/161627-horn-honk-wanted.html
%
% Diy_Audio_WTF (Wavelet Audio Analysis) 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, or (at your option) any
% later version.
%
% Diy_Audio_WTF 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 Diy_WAA; see the file COPYING. If not, write to the Free
% Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
%
%--------------------------------------------------------

if d_bug == 1;
disp ("DEBUG REPORT: processing function : wtf_mag_norm , please WAIT !" )
endif

max_wtf_mag = max(max(wtf_mag));
min_wtf_mag = min(min(wtf_mag));

wtf_mag_norm = wtf_mag - max_wtf_mag; % set matrix peak value to 0

if d_bug == 1;
disp ("DEBUG REPORT: function : wtf_mag_norm processed successfully" )
endif

endfunction

mige0 · 2010-04-19 9:58 pm

function [wtf_t, wtf_f, wtf_mag] = wtf_lecleach (ir_time_trim, ir_amp_trim, d_bug);
%
% This file is part of Diy_Audio_WTF package, intended to get started on Audio Wavelet Transformation
% on free Octave software
% by Jean-Michel Le Clac'h (version 2010)
% adapted by Michael Gerstgrasser 2010
%
% "wtf_lecleach" is a function to perform Jean-Michel LeCleach spectrogram wavelt analysis.
% "wtf_lecleach" is EXCLUDED from GNU licence of the Diy_Audio_WTF package.
% Ask original Author for permission of any further usage.
%
% Input variables :
% ir_time_trim is the vector containing the time data of the IR file
% ir_amp_trim is the vector containing the amplitude of the IR file
% d_bug is for providing additional information in case of debug mode
%
% Output variables:
% wtf_x time scale variable
% wtf_y frequency scale variabel
% wtf_mag magnitude variable
%
% see also load_wav_ir.m, load_ascii_ir.m, load_ascii1_ir.m, fft_ir_raw, waveletdemo.m

% Contributing authors to the Diy_Audio_WTF package in order of code adaption date:
%
% Michael Gerstgrasser - Contact information: user 'mige0' at www.diyaudio.com
% Elias Pekonen - Contact information: user 'Elias' at www.diyaudio.com
% J. Lewalle - Contact information: Prof. Jacques Lewalle
% Jean-Michel Le Clac'h - Contact information: user 'Jmmlc' at www.diyaudio.com
%
% check out on www.diyaudio.com :
% "Help on making MULTIRESOLUTION WAVELET ANALYSIS available ?"
% http://www.diyaudio.com/forums/soft...ltiresolution-wavelet-analysis-available.html
% "WTF!? Wavelet TransForm for audio measurements - What-is? and How-to?"
% http://www.diyaudio.com/forums/mult...et-transform-audio-measurements-what-how.html
% "Horn Honk $$ WANTED $$"
% http://www.diyaudio.com/forums/multi-way/161627-horn-honk-wanted.html
%
% Diy_Audio_WTF (Wavelet Audio Analysis) 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, or (at your option) any
% later version.
%
% Diy_Audio_WTF 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 Diy_WAA; see the file COPYING. If not, write to the Free
% Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
%
%--------------------------------------------------------

if d_bug == 1;
disp ("DEBUG report: processing function : wtf_lecleach , please WAIT !" )
endif

%
%=================================================================
% adapt interfacing of variables
%=================================================================

sig=ir_amp_trim;
N=length(sig);
s_r_rec=10/(ir_time_trim(110)-ir_time_trim(100)); % reconstruckt sample rate
Fe=s_r_rec;
Fy=Fe;
period=1/Fe;
%
%=================================================================
% SIGNAL PRECONDITIONNING
%=================================================================
% limits and number of samples of the signal window
%=================================================================
%
fraction=0.3;
% 12 periods ==> fraction = 0.17
tc = gauspuls2('cutoff',600.0,fraction,[],-40);
tc=period*round(tc/period);
tp = -tc : period : tc;
Npts=length(tp);
y=zeros(1,Npts);
%
%=================================================================
% centering of the Impulse Response inside the signal window
%=================================================================
%
ymax=max(sig); ymin=min(sig);
imax= find(sig==ymax);
if ymax<abs(ymin);
imax=find(sig==ymin);
endif
idec=round(Npts-1)/2;
idec= idec-imax+1;
for i = 1:N;
j=i+idec;
if j>0;
if j<=Npts;
y(j)=sig(i);
endif
endif
endfor
%===============================================================
% gaussian envelop pulse with frequency center = 1000Hz
%=================================================================
%
fcentr=1000;
tc = gauspuls2('cutoff',fcentr,2.0,[],-40); % tc = cut off time
tc=period*round(tc/period);
t = -tc : period : tc;
[cosgauss, singauss, envelop] = gauspuls2(t, fcentr, 2.0);
%
%===============================================================
% convolution of the IR by the gaussian envelop pulse
%===============================================================
%
zcos=conv(y,cosgauss);
zsin=conv(y,singauss);
zcos2 = zcos.*zcos;
zsin2 = zsin.*zsin;
zamp=zcos2+zsin2;
zamp = zamp.^(1/2);
long_decal= round((length(t)+0.5)/2);
ymax=max(zamp);
if ymax>0;
imax= find(zamp==ymax);
else
imax=round((Npts+0.5)/2);
endif
imax=imax-long_decal+1-round((Npts+0.5)/2);

%===============================================================
% limits of frequency of the frequency analysis
%===============================================================

fmax=23000;
if fmax>Fe/2;
fmax=Fe/2;
endif
Nfreqs=400;
fmin=200;
if fmin<10;
fmin=10;
endif
freq=ones(1,Nfreqs);
sfen=ones(1,Nfreqs);
logfreq=zeros(1,Nfreqs);
freq(1)=fmin;
dfreq=(fmax/fmin)^(1/(Nfreqs-1));
for i = 2:Nfreqs
freq(i)=freq(i-1)*dfreq;
endfor
logfreq=log10(freq); % for vertical axis of the spectrogram

%===============================================================
%===============================================================
% MAIN LOOP
%===============================================================
%
spectro=zeros(Nfreqs,Npts);

for i = 1:Nfreqs
%===============================================================
% calculation of a gaussian envelop pulse for the ith frequency
%===============================================================
fcentr= freq(i);
tc = gauspuls2('cutoff',fcentr,1.0,[],-40);
tc=period*round(tc/period);
t = -tc : period : tc;
[cosgauss,singauss,envelop] = gauspuls2(t,fcentr,1.0);
%===============================================================
% convolution of the IR by the gaussian envelop pulse
%===============================================================
zcos=conv(y,cosgauss); zsin=conv(y,singauss);
zcos2 = zcos.*zcos;
zsin2 = zsin.*zsin;
zamp=zcos2+zsin2;
zamp = zamp.^(1/2);
long_decal= round((length(t)+0.5)/2);
spectro(i,

=zamp(long_decal:long_decal+Npts-1);
endfor

%===============================================================
% spectrogram in dBs
%===============================================================

dBmin=-40; % 40 dB color scale
spectro = 20*log10(spectro);
spectro=max(dBmin,spectro);

%===============================================================
% graphical output of the spectrogram
%===============================================================

spectro1=spectro

, 1);
f = (1:length (spectro1));
tp=tp*1000;

%===============================================================
% output for further processing
%===============================================================

wtf_t = tp;
wtf_f = freq;
wtf_mag = spectro;

if d_bug == 1;
disp ("DEBUG report: function : wtf_lecleach processed successfully" )
endif

%===============================================================
% END
%===============================================================
endfunction

mige0 · 2010-04-19 9:59 pm

function [wtf_t, wtf_f, wtf_mag] = wtf_morlet_lewalle (ir_time_trim, ir_amp_trim, d_bug);

%--------------------------------------------------------
% wtf_morlet_lewalle
%--------------------------------------------------------
% This file is part of Diy_Audio_WTF package, intended to get started on Audio Wavelet Transformation
% on free Octave software
% Copyright (C) J. Lewalle 1995
% rewritten by Michael Gerstgrasser 2010
%
% "wtf_morlet_lewalle" is a function for continuous wavelet processing
%
% Input variables :
% ir_time_trim is the vector containing the time data of the IR file
% ir_amp_trim is the vector containing the amplitude of the IR file
% d_bug is for providing additional information in case of debug mode
%
% Output variables:
% wtf_x time scale variable
% wtf_y frequency scale variabel
% wtf_mag magnitude variable
%
%
% see also load_wav_ir.m, load_ascii_ir.m, load_ascii1_ir.m, fft_ir_raw, waveletdemo.m

% Contributing authors to the Diy_Audio_WTF package in order of code publishing date:
%
% Michael Gerstgrasser - Contact information: user 'mige0' at www.diyaudio.com
% Elias Pekonen - Contact information: user 'Elias' at www.diyaudio.com
% J. Lewalle
%
% check out on www.diyaudio.com :
% "Help on making MULTIRESOLUTION WAVELET ANALYSIS available ?"
% http://www.diyaudio.com/forums/soft...ltiresolution-wavelet-analysis-available.html
% "WTF!? Wavelet TransForm for audio measurements - What-is? and How-to?"
% http://www.diyaudio.com/forums/mult...et-transform-audio-measurements-what-how.html
% "Horn Honk $$ WANTED $$"
% http://www.diyaudio.com/forums/multi-way/161627-horn-honk-wanted.html
%
% Diy_Audio_WTF (Wavelet Audio Analysis) 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, or (at your option) any
% later version.
%
% Diy_Audio_WTF 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 Diy_WAA; see the file COPYING. If not, write to the Free
% Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
%

if d_bug == 1;
disp ("DEBUG REPORT: processing function : wtf_morlet_lewalle , please WAIT !" )
endif

%--------------------- interface with input variables -------------------------------

dx = (ir_time_trim(110)-ir_time_trim(100)) / 10;

nvcs = 20 ; % number of voices per octave
drs = 2^(2/nvcs);
nx = length(ir_amp_trim);
ns = log2(nx)*nvcs;
% covers the range of time-scales in the signal;
% can be narrowed to the range of interest.
ns = (fix(ns/2))*2+1; % make sure ns is odd integer

%--------------------- removes the mean from a signal -------------------------------

y = ir_amp_trim - mean(ir_amp_trim);

%-------------------- START script ----------------------------------------------------------------

nx = length

;
x=linspace(0,nx-1,nx)*dx;
ff = lewa_fftfreq(x,dx, d_bug); % Fourier frequency axis (2-sided linear)
yh = fft

; % FFT of signal

s(ns) = (dx/pi)^2; % dx/pi is Nyquist cutoff
% see comment regarding energy in the `showmaps' function
% see the extension of the spectra on the figure.
for is=ns-1:-1:1;
s(is) = s(is+1)*drs;
endfor
for is=1:ns ;
kf(is) = sqrt(1/s(is))/2/pi;
endfor

% ----------------- MORLET TRANSFORM: norm, phase and complex ---------------------------

wtf_mag = lewa_mymorlet(yh,ff,kf, d_bug);
wtf_mag = log10(wtf_mag) .*20;
wtf_t = x*1000; % scale in ms
wtf_f = kf;

if d_bug == 1;
disp ("DEBUG REPORT: function : wtf_morlet_lewalle processed successfully" )
endif

endfunction

mige0 · 2010-04-19 10:00 pm

function [ff] = lewa_fftfreq(x,dx, d_bug)

%--------------------------------------------------------
% lewa_fftfreq
%--------------------------------------------------------
% This file is part of Diy_Audio_WTF package, intended to get started on Audio Wavelet Transformation
% on free Octave software
% Copyright (C) J. Lewalle 1995
% rewritten by Michael Gerstgrasser 2010
%
% " lewa_fftfreq" constructs the linear 2-sided scale of frequencies for FFT
%
% Input variables :
% d_bug is for providing additional information in case of debug mode
%
%
% see also load_wav_ir.m, load_ascii_ir.m, load_ascii1_ir.m, fft_ir_raw, waveletdemo.m

% Contributing authors to the Diy_Audio_WTF package in order of code publishing date:
%
% Michael Gerstgrasser - Contact information: user 'mige0' at www.diyaudio.com
% Elias Pekonen - Contact information: user 'Elias' at www.diyaudio.com
% J. Lewalle
%
% check out on www.diyaudio.com :
% "Help on making MULTIRESOLUTION WAVELET ANALYSIS available ?"
% http://www.diyaudio.com/forums/soft...ltiresolution-wavelet-analysis-available.html
% "WTF!? Wavelet TransForm for audio measurements - What-is? and How-to?"
% http://www.diyaudio.com/forums/mult...et-transform-audio-measurements-what-how.html
% "Horn Honk $$ WANTED $$"
% http://www.diyaudio.com/forums/multi-way/161627-horn-honk-wanted.html
%
% Diy_Audio_WTF (Wavelet Audio Analysis) 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, or (at your option) any
% later version.
%
% Diy_Audio_WTF 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 Diy_WAA; see the file COPYING. If not, write to the Free
% Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
%

if d_bug == 1;
disp ("DEBUG REPORT: processing function : lewa_fftfreq , please WAIT !" )
endif

n=length(x);
for ii=1

n/2) ; ff(ii)= (ii-1)/n/dx; end
for ii=(1+n/2):n ; ff(ii)= -(n-ii+1)/n/dx; end

if d_bug == 1;
disp ("DEBUG REPORT: function : lewa_fftfreq processed successfully" )
endif

endfunction

mige0 · 2010-04-19 10:00 pm

function [wtf_mag] = lewa_mymorlet(xh,ff,kf, d_bug)
% xh: DFFT of signal
% ff: vector of frequencies (linear)
% kf: vector of frequencies (logarithmic)
% z0: parameter
% yn: norm of Morlet transform arrays
%--------------------------------------------------------
% lewa_mymorlet
%--------------------------------------------------------
% This file is part of Diy_Audio_WTF package, intended to get started on Audio Wavelet Transformation
% on free Octave software
% Copyright (C) J. Lewalle 1995
% rewritten by Michael Gerstgrasser 2010
%
% "lewa_mymorlet" is a function for Morlet wavelet trasnform of a signal
%
% Input variables :
% d_bug is for providing additional information in case of debug mode
%
%
%
%
% see also load_wav_ir.m, load_ascii_ir.m, load_ascii1_ir.m, fft_ir_raw, waveletdemo.m

% Contributing authors to the Diy_Audio_WTF package in order of code publishing date:
%
% Michael Gerstgrasser - Contact information: user 'mige0' at diyAudio
% Elias Pekonen - Contact information: user 'Elias' at diyAudio
% J. Lewalle
%
% check out on diyAudio :
% "Help on making MULTIRESOLUTION WAVELET ANALYSIS available ?"
% http://www.diyaudio.com/forums/soft...ltiresolution-wavelet-analysis-available.html
% "WTF!? Wavelet TransForm for audio measurements - What-is? and How-to?"
% http://www.diyaudio.com/forums/mult...et-transform-audio-measurements-what-how.html
% "Horn Honk $$ WANTED $$"
% http://www.diyaudio.com/forums/multi-way/161627-horn-honk-wanted.html
%
% Diy_Audio_WTF (Wavelet Audio Analysis) 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, or (at your option) any
% later version.
%
% Diy_Audio_WTF 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 Diy_WAA; see the file COPYING. If not, write to the Free
% Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
%
%--------------------------------------------------------

if d_bug == 1;
disp ("DEBUG REPORT: processing function : lewa_mymorlet , please WAIT !" )
endif

n=length(xh);
ns=length(kf);
sq2pi = sqrt(2*pi);
z0 = 5; % variable
cps= 1.4406;
for is=1:ns
for ii=1:n
z=ff(ii)/kf(is);
morh(ii)= z0/kf(is)/sq2pi * (exp(-(z0^2*(1-z)^2)/2) ...
- exp(-z0^2*(1+z^2)/2) ) +0*i;
morh(ii) = morh(ii) * kf(is)/sqrt(cps) ;% my scaling
endfor
for ii=1:n ;
yh(ii) = xh(ii)*morh(ii);
endfor
ycmor(is,

= ifft(yh); % the complex inverse transform
end
wtf_mag = abs(ycmor); % its norm

if d_bug == 1;
disp ("DEBUG REPORT: function : lewa_mymorlet processed successfully" )
endif

endfunction

mige0 · 2010-04-19 10:02 pm

to be continued....

by a zip file for download of complete package

have fun !

Michael

mige0 · 2010-04-21 6:31 am

To install Octave on windows

1.) download and install
Octave-3.2.4_i686-pc-mingw32_gcc-4.4.0_setup.exe

2.) download and install
Xoctave_setup 0.9.92.exe

3.) Start XOctave

4.) on first call Xoctave tells you the Octave path is not assigned
>> view >> options >> console >> "where is Octave"
enter your path to Console software (most probably: "C:\Octave\3.2.4_gcc-4.4.0\bin\octave.exe")

5.) download DIY_wavelet_Package and unpack
place folder "mige0_start_wavelet_ver_xx" to wherever you like

6.) tell Octave the path where your files are to be found
>> view >> path >>
enter your path to "mige0_start_wavelet_ver_xx" folder and press "Add Path"
check if your path to "mige0_start_wavelet_ver_xx" has been assigned by pressing "Update"

now >> view >> Console >> and you are ready to go under Windows !

################

To install Octave under Ubuntu (10.04 beta I've tryed), its already listed under Software and you can install comfortably both, Octave and QtOctave from there.

You do not have to tell QtOctave where to find Octave but you have to tell QtOctave where to find "mige0_start_wavelet_ver_xx" folder of course.

Finally set cursor to " console " and you are ready to go under Linux !

################

type in " start_wavelet " and follow the user menu you are told – or – check out Elias' " waveletdemo "

For convenience I have added an IR_file.wav file from my 0-deg honker – so you can go through the menu to check out easily.

If you are going to analyse your own IR files, you either have to place it in the "mige0_start_wavelet_ver_xx" folder, or you have to add an additional path to your preferred directory.

################

Navigating through the Diy_Audio_WTF (Wavelet Audio Analysis) package:

>>> start_wavelet

will start the WTF package and will bring up the intro and in the first menu you will be ask to choose between :

- normal operation
- debug operation

I recommend to use "normal operation" unless you face unexpected problems.
"debug operation" gives you additional information on what functions are called and have processed successfully. So this is a handy tool to post a simple screenshot in case something went wrong.

********

in next menu you will be asked to select an impulse response file.

You can just press enter to select the default file "ir_file.wav" of type in the name of your file.
Its best to avoid spaces and any special characters for your impulse response file.
Also the file has to be located where it can be found by octave – meaning – you either can put it into the "mige0_start_wavelet_ver_xx" folder, or you set an additional path to your preferred folder.

********

in next menu you will be asked to select an impulse response file type.

As for now the package can process different file types:
- WAV file
- ARTA text file export
- plain two column ascii text file
- plain one column ascii text file

in case of "ascii two column text file" it is assumed that the first column contains the time and the second column contains the amplitude of your impulse response.

in case of "ascii one column text file" it is assumed that the first and only column contains the amplitude of your impulse response and you will be asked to enter the sample rate in a sub menue.

********

in next menu you will be asked to select wavelet analysis type.

As for now the package can process the "Lewalle Morlet wavelet analysis" only, but is fully prepared for the Jean-Michel Le Clea'ch wavelet analysis" as well.

********

The first plot popping up after the third menu will display the raw impulse response from 0 to 100ms

The second plot popping up after the fourth menu will display the normalized impulse
response from 0 to 13ms with the peak centred at 3ms

The third plot popping up after the fourth menu will display the wavelet analysis for a range of 40dB and within a time range of 0 to 13ms with the peak centred at 3ms

The fourth plot popping up after the fourth menu will display the normalized wavelet analysis for a range of 40dB and within a time range of 0 to 13ms with the peak centred at 3ms

********

Between the menus additional information is displayed that you may you find interesting and helpful.

If you quit by " EXIT " in menu, you have to reconnect Xoctave to Octave.
Be warned that for one reason or another, the plots are pretty slow when you work under Windows.

********

Hope you have fun with the first free software package available, aimed for wavelet audio analysis to the "non-math-head" user.

You can download the full package here:

http://www.kinotechnik.edis.at/pages/diyaudio/wavelet/mige0_start_wavelet_ver_1_0.zip

In the "readme" you find above text, and I also included the screen shots of the menu structure and the plots that get calculated.

Michael

maiky76 · 2010-04-21 7:49 am

Hi,

Thanks for the contribution!
In case you don't know there is this amazing Octave UI:
QtOctave: a front-end for Octave/un Front-End para Octave
I have been using it for the past year and like it better than Xoctave.
It is compatible with many Linux flavors and Windows.

Cheers
Mickael

maiky76 · 2010-04-21 8:26 am

By the way whay using a none free app such as Arta when a package is especialy designed for Matlab/Octave.
MATAA (Mat's Audio Analyzer) | Get MATAA (Mat's Audio Analyzer) at SourceForge.net
I also use it and it is fantastic and since you guys are familiar with Matlab/Octave and signal processing it should be ok.
I sucessfully used it with an Emu 0404 USB and a M-audio Transit, I coudl help to get you started if you want.

Gaga · 2010-04-21 10:32 am

Wavelet Analysis for OS X Platform

Dear all,

Though I do not understand everything on wavelet-analysis discussed here, I followed the threads on 'how to' and horn honk with great interest.

The main contributors here seem to work on windows or linux platforms. As an os x user I just want to direct you to Christoph Lauers work and homepage on audio and wavelet analysis (www.christoph-lauer.de) - may be you missed his nice work for this reason.

Please refer to 'Sonogram' and the 'Fast-Wavelet-Transformation Toolbox in C++' (see the Blog on his homepage). May be of interest anyway.

Regards,
gaga

doug20 · 2010-04-21 3:05 pm

It works, Thank you so much!!

mige0 · 2010-04-21 6:26 pm

doug20 said:
It works, Thank you so much!!

Good !

In case you haven't already corrected for yourself, I fixed the minor bug in the plotting x-axis of the wavelets.

You can download the zip file again and replace "start_wavelet.m" if you like (didn't do an new version for that)

Michael

mige0 · 2010-04-21 6:31 pm

maiky76 said:
By the way whay using a none free app such as Arta when a package is especialy designed for Matlab/Octave.
MATAA (Mat's Audio Analyzer) | Get MATAA (Mat's Audio Analyzer) at SourceForge.net
I also use it and it is fantastic and since you guys are familiar with Matlab/Octave and signal processing it should be ok.
I sucessfully used it with an Emu 0404 USB and a M-audio Transit, I coudl help to get you started if you want.

Thanks for pointing there !

maiky76 said:
Hi,

Thanks for the contribution!
In case you don't know there is this amazing Octave UI:
QtOctave: a front-end for Octave/un Front-End para Octave
I have been using it for the past year and like it better than Xoctave.
It is compatible with many Linux flavors and Windows.

Cheers
Mickael

I missed that thanks, its good to have QtOctave under windows – way better than XOctave I certainly agree.
Those Spanish speaking guys seem to be pretty active. Lets hope they do a "monkey interface" one time as well....

mige0 · 2010-04-21 6:41 pm

Gaga said:
Dear all,

Though I do not understand everything on wavelet-analysis discussed here, I followed the threads on 'how to' and horn honk with great interest.

The main contributors here seem to work on windows or linux platforms. As an os x user I just want to direct you to Christoph Lauers work and homepage on audio and wavelet analysis (www.christoph-lauer.de) - may be you missed his nice work for this reason.

Please refer to 'Sonogram' and the 'Fast-Wavelet-Transformation Toolbox in C++' (see the Blog on his homepage). May be of interest anyway.

Regards,
gaga

Looks very interesting indeed.

But too much troubles to get involved into C++ coding - actually this "DIY_WTF_package" is only a detour project of mine, to have a tool available to refine development of Gaussian MinPhase dipole horns - already quite a huge detour actually (in terms of "quality time consumption") - though it was a lot of fun too...

Its been a while that I did some more serious coding in PASCAL - LOL - I guess it must be roughly two decades now...

Michael

doug20 · 2010-04-21 6:55 pm

mige0 said:
Good !

In case you haven't already corrected for yourself, I fixed the minor bug in the plotting x-axis of the wavelets.

You can download the zip file again and replace "start_wavelet.m" if you like (didn't do an new version for that)

Michael

I deleted the original code with the latest zip, that should correct it, right?

mige0 · 2010-04-21 7:24 pm

Anybody who can help in refining ?

Especially the color bar isn't very pretty compared to the slim presentation in CLIO :

http://www.diyaudio.com/forums/multi-way/161627-horn-honk-wanted-12.html#post2154127

Also the grid of my plots could need some brush up...

Michael

Elias · 2010-04-22 12:26 pm

Impulse response file conversion CLIO .mls -> WAV -> Octave

Hello,

Here I attach an Octave script to convert CLIO exported impulse response file to WAV file. Octave can read WAV file directly.

Note! This conversion method is purely experimental, since I have no information on CLIO .mls binary file header syntax. However, it works if you enter parameters manually.

The parameters are:
* impulse response file name and path
* sampling frequency
* start and stop indexes

If the file header length remains the same Nstart = 250 is allways valid. However, Nstop depends on impulse response length and is known only by the original user, and thus one may need to change Nstop value.

One plot of the imported impulse response is generated in dB scale to determine if import is succesfull.

I have tested it with ZilchLab's impulse response file from here:
http://www.diyaudio.com/forums/multi-way/161627-horn-honk-wanted-13.html#post2155094

Feel free to attach this to the diyaudio wavelet package

Code:

### Reads exported binary impulse response file
### and converts it to WAV file
### and plots the impulse response
###
### file version 1.1
###
### (c) 2010 Elias Pekonen  
###
### 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 3 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 <http://www.gnu.org/licenses/>.


###########################################################
### ENTER PARAMETRS HERE

### NAME OF THE IMPULSE RESPONSE FILE
filename = 'C:/your_folder_path/your_file.mls'


### SAMPLING FREQUENCY
Fs = 48000


### PARAMETERS TO SKIP BINARY FILE HEADER AND TAIL
Nstart = 250
Nstop = 16000



###########################################################
fid = fopen (filename);

[val, count] = fread (fid, Inf, 'float32');

fclose (fid);

ir = val (Nstart : min (Nstop, count));

wavwrite (ir, Fs, strrep(filename, '.mls', '.wav'));

figure (1)
plot ((0 : length(ir) - 1)/Fs, 20*log10 (abs (ir)));
xlabel ('Time [s]');

- Elias

Elias · 2010-04-22 12:53 pm

Hi,

It may have something to do which option you have selected as your graphic backend. There are different versions of the GNU Plot for example. Octave uses GNU Plot to display images

- Elias

mige0 said:
Anybody who can help in refining ?

Especially the color bar isn't very pretty compared to the slim presentation in CLIO :
...

Also the grid of my plots could need some brush up...

Michael

mige0 · 2010-04-22 9:03 pm

doug20 said:
I deleted the original code with the latest zip, that should correct it, right?

Yes - should work...

###########

Elias said:
Here I attach an Octave script to convert CLIO exported impulse response file to WAV file.
...

Feel free to attach this to the diyaudio wavelet package

- Elias

Thanks. Elias !
Works fine !

in case anybody would like to analyse CLIO files, just check out these intermediate versions below:

Code:

%--------------------------------------------------------
%	start_wavelet
%--------------------------------------------------------	
%	This file is part of Diy_Audio_WTF package, intended to get started on Audio Wavelet Transformation
%	on free Octave software  
%	Copyright (C) 2010  Michael Gerstgrasser
% 	
%	"start_wavelet" is a script to load impulse response files of various formats
%	and to check valid import by simple impulse response and frequency response plot 
%
%	ir_time is the vector containing the time dat	a of the IR file
%	ir_amp is the vector containing the amplitude of the IR file
%	
%
%	see also :
%	start_wavelet.m, load_arta_txt_ir.m, load_ascii_ir.m, load_ascii1_ir.m, 
%	load_wav_ir.m, shape_ir.m, matrix_norm.m, wtf_mag_norm.m, wtf_lecleach.m, 
%	wtf_morlet_lewalle.m,  lewa_fftfreq.m, lewa_mymorlet.m, waveletdemo.m, content.m,
%
%	Contributing authors to the Diy_Audio_WTF package in order of code publishing date:
%
%	Michael Gerstgrasser 	- Contact information:  user 'mige0' at www.diyaudio.com
%	Elias Pekonen 		- Contact information:  user 'Elias' at www.diyaudio.com
%	Jacques Lewalle  	- http://www.ecs.syr.edu/faculty/lewalle/
%	Jean-Michel Le Clac'h	- Contact information:  user 'Jmmlc' at www.diyaudio.com
%
%
% 	Diy_Audio_WTF  (Wavelet Audio Analysis) package 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, or (at your option) any      
% 	later version. 
%
%                                                    
% 	check out on www.diyaudio.com :
%	"Help on making MULTIRESOLUTION WAVELET ANALYSIS available ?"
%	http://www.diyaudio.com/forums/software-tools/163473-help-making-multiresolution-wavelet-analysis-available.html
%	"WTF!? Wavelet TransForm for audio measurements - What-is? and How-to?"
%	http://www.diyaudio.com/forums/multi-way/164029-wtf-wavelet-transform-audio-measurements-what-how.html
%	"Horn Honk $$ WANTED $$"
%	http://www.diyaudio.com/forums/multi-way/161627-horn-honk-wanted.html
%                                                          
%                                                                          
% 	Diy_Audio_WTF 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 Diy_WAA; see the file COPYING.  If not, write to the Free    
% 	Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.             
%                                                                          
%--------------------------------------------------------

%     
%--------------------------------------------------------
% 	global declaration    
%--------------------------------------------------------

global ir_file_name;

%--------------------------------------------------------
% 	introduction and prompt for debug option
%--------------------------------------------------------
disp("")
disp("")
disp("#######+++++++++++++++++++++++++++++++++++++++++##########")
disp("")
disp("")
disp("HELLO !     Wecome to DIY Audio Wavelet Analysis ! ")
disp("contributing authors : ")
disp("Michael Gerstgrasser ")
disp("Elias Pekonen")
disp("Jaques Lewalle ")
disp("Jean-Michel Le Cleac'h")
disp("")
disp("You will be promted to pick your choices... ")
disp("")

k=0;
d_bug=1;
k=menu("Enter 1 for normal operation"," normal operation"," debug mode", " EXIT");

if k==1;
	d_bug = 0;
endif;
if k==2;
	d_bug=1;
	disp("")
	disp("DEBUG mode selected !  Additional progress information will be displayed")
endif;
if k==3;
	disp("")
	disp("DIY wavlet analysis will be terminated")
	quit;
endif;	


%--------------------------------------------------------
% 	prompt for file name
%--------------------------------------------------------

disp("")
disp("Make sure your impulse response file is located where it can be found by Octave")
disp("Make sure your impulse response file name avoids spaces and special characters")
disp("")
fflush (stdout);
ir_file_name=input("HELLO !  Please type in file name ( press ENTER for ir_file.wav )","s");

	if isempty(ir_file_name);
	ir_file_name="ir_file.wav";			% assign default
	endif;
	
ir_file_name						% display selected file name
disp("")
disp("")



%--------------------------------------------------------
% 	prompt for file type
%--------------------------------------------------------

k=0;


k=menu("HELLO !  Please select your impulse response file type","32 bit or 16 bit WAV file","ARTA txt file","two column TEXT file","one column TEXT file","CLIO import ( experimental )","EXIT");



	if k==1,
		disp("")	
		disp("please wait... ")
		disp(ir_file_name)						% display selected file name
		disp("16 bit or 32 bit WAV file will be imported")
		disp("")	
		disp("Make sure to use a 16 bit or 32 bit WAV impulse respones file")
		[ir_time ir_amp] = load_wav_ir(d_bug);		% load wav32 impulse response file and assign to variables ir_time, ir_amp
	endif
	if k==2,
		disp("")	
		disp("please wait... ")
		disp(ir_file_name)	
		disp("ARTA text file  will be imported")
		disp("")	
		disp("Make sure to use a TEXT impulse respones file exported from ARTA")
		[ir_time ir_amp] = load_arta_txt_ir(d_bug);	% load ARTA txt impulse response file and assign to variables ir_time, ir_amp
	endif
	if k==3, 
		disp("")	
		disp("please wait... ")
		disp(ir_file_name)						% display selected file name
		disp("Two column TEXT file will be imported")
		disp("")	
		disp("Make sure to use a plain TEXT, headerless, two column impulse response file")
		disp("Make sure first column in your impulse response file contains time")
		disp("Make sure second column in your impulse response file contains amplitude")
		[ir_time ir_amp] = load_ascii_ir(d_bug);			% load ascii impulse response file and assign to variables ir_time, ir_amp
	endif
	if k==4, 
		disp("")
		disp("please wait... ")	
		disp(ir_file_name)						% display selected file name
		disp("One column TEXT file will be imported")
		disp("")	
		disp("Make sure to use a plain TEXT, headerless, one column impulse response file")
		disp("Make sure first column in your impulse response file contains amplitude")
		disp("Check the following plot for correct impulse response file import")
		[ir_time ir_amp] = load_ascii1_ir(d_bug);			% load ascii impulse response file and assign to variables ir_time, ir_amp
	endif
	if k==5, 
		disp("")
		disp("please wait... ")	
		disp(ir_file_name)						% display selected file name
		disp("One column TEXT file will be imported")
		disp("")	
		disp("Make sure to use a plain TEXT, headerless, one column impulse response file")
		disp("Make sure first column in your impulse response file contains amplitude")
		disp("Check the following plot for correct impulse response file import")
		[ir_time ir_amp] = load_clio_ir(d_bug);			% load ascii impulse response file and assign to variables ir_time, ir_amp
	endif
	if k==6,
		disp("")
		disp("DIY wavlet analysis will be terminated")
		quit;
	endif;
%--------------------------------------------------------
%	restore  and display : samplerate , sample count , total time
%--------------------------------------------------------

disp("")
disp("")
disp("HELLO ! Properties of imported file :")
disp("")
sample_rate = 10/(ir_time(110)-ir_time(100))			% reconstruckt sample rate
samples = length(ir_amp)							% sample count			
total_time = samples / sample_rate					% calculate total time of measurement


%--------------------------------------------------------
%--------------------------------------------------------


%--------------------------------------------------------
%	plot impulse response to check correct import
%--------------------------------------------------------
clf;
disp("")
disp("")
disp("HELLO !  Please check the following plot for correct impulse respones file import")


figure(1)
plot (ir_time,ir_amp,"-;Impulse Response;")
set (gca (), "xlim", [0 0.1]);
title(["HELLO !   Please check plot below for correct impulse response import (first 100ms) ! "," Imported file: ",(ir_file_name)])

xlabel ("Time in seconds");
ylabel ("Amplitude");



samples_trim = 1024;	
ir_time_trim=ir_time;
ir_amp_trim=ir_amp;
[ir_time_trim, ir_amp_trim] = shape_ir (ir_time, ir_amp, samples_trim, d_bug);


%--------------------------------------------------------
%	select wavelet analysis
%--------------------------------------------------------



k=0;
k=menu("HELLO ! please select wavlet analysis"," Jean-Michel Le Cleac'h Wavelet Analysis"," Jaques Lewalle Morlet Wavelet Analysis", " EXIT");

if k==1;
	disp("")
	disp("HELLO !  Wavelet analysis in progress - this may take a few minutes")
	disp("------------	normalised impulse response will be displayed	------------  ")
	disp("------------		wavelet plot will be displayed		------------  ")
	disp("------------	normalised wavelet plot will be displayed		------------  ")
	[wtf_t, wtf_f, wtf_mag] = wtf_lecleach (ir_time_trim, ir_amp_trim, d_bug) ;
endif;
if k==2;
	disp("")
	disp("HELLO !  Wavelet analysis in process - this may take a few minutes")
	disp("------------	normalised impulse response plot will be displayed		------------  ")
	disp("------------		wavelet plot will be displayed		------------  ")
	disp("------------	normalised wavelet plot will be displayed		------------  ")
	[wtf_t, wtf_f, wtf_mag] = wtf_morlet_lewalle (ir_time_trim, ir_amp_trim, d_bug);
endif;
if k==3;
	disp("")
	disp("DIY wavelet analysis will be terminated")
	quit;
endif;	


%--------------------------------------------------------
%	plot normalised impulse response 
%--------------------------------------------------------




figure(2)
clf(2);
plot (ir_time_trim,ir_amp_trim,"-;Impulse Response;")
set (gca (), "xlim", [0 0.1]);
title(["HELLO !   Please check plot below for correct impulse response import (first 100ms) ! "," Imported file: ",(ir_file_name)])
set (gca (), "xlim", [0  0.013]);
xlabel ("Time in seconds");
ylabel ("Amplitude");


%--------------------------------------------------------
%	plot  wavlet analysis
%--------------------------------------------------------

fmin = 1000;											%	or check for min freqeuncy to plot	min(wtf_f); 
fmax = max (wtf_f);									%	check for max freqeuncy to plot
mag_low_limit=-40;									%	set minimum magnitude


wtf_t = wtf_t - 3;									% 	set max peak in wavelet plot to zero time
max_wtf_mag = max(max(wtf_mag));						% 	set peak magnitude to 0 dB
wtf_mag_norm = wtf_mag - max_wtf_mag;		
wtf_mag_norm = max(mag_low_limit,wtf_mag_norm);		% 	set minimum magnitude to mag_low_limit


figure (3);
clf(3);
pcolor (wtf_t, wtf_f, wtf_mag_norm);		
shading ("interp");
 colorbar ();											%	optional : colorbar ();
xlabel ("Time [ms]");
ylabel ("Freq [Hz]");
set (gca (), "yscale", "log");
set (gca (), "xlim", [-3 10]);
set (gca (), "ylim", [fmin fmax]);
title(["WTF spectrogram    ",(ir_file_name)])
grid("on");



%--------------------------------------------------------
%	plot normalised wavlet analysis
%--------------------------------------------------------


matrixx_norm = matrix_norm(wtf_mag, d_bug);					%	normalization feature

 
max_wtf_mag = max(max(matrixx_norm));						% 	set peak magnitude to 0 dB
matrixx_norm = matrixx_norm - max_wtf_mag;	 
matrixx_norm = max(mag_low_limit,matrixx_norm);				% 	set minimum magnitude to mag_low_limit
 
 
figure (4);
clf(4);								
pcolor (wtf_t, wtf_f, matrixx_norm);		
shading ("interp");
colorbar ();												% optional : colorbar ();
xlabel ("Time [ms]");
ylabel ("Freq [Hz]");
set (gca (), "yscale", "log");
set (gca (), "xlim", [-3 10]);
set (gca (), "ylim", [fmin fmax]);
title(["WTF spectrogram normalized   ",(ir_file_name)])
grid("on");

http://www.kinotechnik.edis.at/page...e0_start_wavelet_ver_inter1_1/start_wavelet.m

Code:

function [ir_time, ir_amp] = load_clio_ir (d_bug);	### Reads exported binary impulse response file
### and converts it to WAV file
### and plots the impulse response
###
### file version 1.1
###
### (c) 2010 Elias Pekonen  
###
### 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 3 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 <http://www.gnu.org/licenses/>.


###########################################################
### ENTER PARAMETRS HERE

### 	%NAME OF THE IMPULSE RESPONSE FILE


%--------------------------------------------------------

if d_bug == 1;
	disp ("DEBUG report: processing function : load_wav_ir , please WAIT !" )
endif
%--------------------------------------------------------
% 	global declaration    
%--------------------------------------------------------

global ir_file_name;

%--------------------------------------------------


%--------------------------------------------------------
% 	prompt for sampling rate
%--------------------------------------------------------

disp("")
disp("NOTE !  This conversion method is purely experimental. Please report back your findings")
disp("")
fflush (stdout);
sampling_freq	= input("HELLO !  Type in sample rate [Hz], or press enter for 48000 Hz");

	if isempty(sampling_freq);
	sampling_freq	= 48000;				% assign default
	endif;

%--------------------------------------------------------

%--------------------------------------------------------
% 	prompt for sampling rate
%--------------------------------------------------------

disp("")
disp("NOTE !  This conversion method is purely experimental. Please report back your findings")
disp("")
fflush (stdout);
sampling_num = input("HELLO !  Enter samle number, or press enter for 16000 samples");

	if isempty(sampling_num);
	sampling_num	= 16000;				% assign default
	endif;

%--------------------------------------------------------


filename = ir_file_name;
ir_file_name 

### SAMPLING FREQUENCY
Fs = sampling_freq
sampling_freq							% display sampling frequency


### PARAMETERS TO SKIP BINARY FILE HEADER AND TAIL
Nstart = 250							% display Nstart
Nstop = sampling_num					% display number of samples



###########################################################
fid = fopen (filename);

[val, count] = fread (fid, Inf, 'float32');

fclose (fid);

ir = val (Nstart : min (Nstop, count));

wavwrite (ir, Fs, strrep(filename, ".mls", ".wav"));
ir_file_name = strrep(filename, ".mls", ".wav")

figure (1)
plot ((0 : length(ir) - 1)/Fs, 20*log10 (abs (ir)));
xlabel ('Time [s]');

[ir_time ir_amp] = load_wav_ir(d_bug);		% load wav32 impulse response file and assign to variables ir_time, ir_amp


endfunction

http://www.kinotechnik.edis.at/page...ge0_start_wavelet_ver_inter1_1/load_clio_ir.m

here is what I got with Zilch's CLIO measurements mentioned above:

#############

Elias said:
Hi,

It may have something to do which option you have selected as your graphic backend. There are different versions of the GNU Plot for example. Octave uses GNU Plot to display images

- Elias

OK

Which one do you use to make the ticks visible on the y-axis.
I didn't get it - also not when I left some spare place - as you did - to not get them overlaid by the plot...

Michael

Elias · 2010-04-23 11:36 am

mige0 said:
Which one do you use to make the ticks visible on the y-axis.
I didn't get it - also not when I left some spare place - as you did - to not get them overlaid by the plot...

Michael

Hi,

I'm using Octave 3.0.3 and GnuPlot 4.3.

How does your pcolor plot looks like with your version if you run the waveletdemo.m?

An externally hosted image should be here but it was not working when we last tested it.

- Elias

Help on making MULTIRESOLUTION WAVELET ANALYSIS available ?

Member

Member

Member

Member

Member

Member

Member

Member

Member

Member

Member

Member

Member

Member

Member

Member

Member

Member

Member

Member