rmeddis@38: function out = mellin_trafo(inx,iny)
rmeddis@38: 
rmeddis@38: %This function computes the Mellin Transformation of a one-dimensional
rmeddis@38: %Signal in analytical terms given as
rmeddis@38: %
rmeddis@38: % S(p)=Integral from 0 to infty of s(t) t^(p-1)dt
rmeddis@38: % Which equals 
rmeddis@38: % S(c) = Integral from -Infty to Infty s(t)*exp(j*c*log(t))d(log(t))
rmeddis@38: %
rmeddis@38: % taken from Irino and Patterson, Speech Communication 2002 Eq. (1) and (3)
rmeddis@38: % Tim Juergens, September 2011
rmeddis@38: %
rmeddis@38: 
rmeddis@38: %Resample the Signal onto a log(t) axis
rmeddis@38: minimalx=min(inx);
rmeddis@38: maximalx=max(inx);
rmeddis@38: logarithmicx= exp([log(minimalx):(log(maximalx)-log(minimalx))/length(inx):log(maximalx)]);
rmeddis@38: logarithmicy= interp1(inx,iny,logarithmicx,'linear','extrap');
rmeddis@38: %figure, semilogx(logarithmicx,logarithmicy);
rmeddis@38: 
rmeddis@38: %Absolute of the inverse Fourier-Transform of the resampled signal
rmeddis@38: out = abs(ifft(logarithmicy));
rmeddis@38: %figure, plot(out(1:40));
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %%just to show that the mellin transform results in invariable patterns if
rmeddis@38: %%the formants are a constant ratio
rmeddis@38: % frequencies_short = [1:8:8000];
rmeddis@38: % frequencies_long = [1:10:10000];
rmeddis@38: % Intervals_long = 1./frequencies_long;
rmeddis@38: % Intervals_short = 1./frequencies_short;
rmeddis@38: % contoursin = sin(2*pi*0.00015.*frequencies_long).^2
rmeddis@38: % figure, plot(frequencies_long,contoursin), hold on, plot(frequencies_short,contoursin,'r')
rmeddis@38: % xlabel('Frequency (Hz)')
rmeddis@38: % ylabel('Rate (arb. units)')
rmeddis@38: % figure, semilogx(Intervals_long,contoursin), hold on, plot(Intervals_short,contoursin,'r')
rmeddis@38: % xlabel('Interval (s)')
rmeddis@38: % ylabel('Rate (arb. units)')
rmeddis@38: % m1 = mellin_trafo(Intervals_long,contoursin);
rmeddis@38: % m2 = mellin_trafo(Intervals_short,contoursin);
rmeddis@38: % figure, plot(m1(1:40)), hold on, plot(m2(1:40),'r');
rmeddis@38: % xlabel('Mellin variable c');
rmeddis@38: % ylabel('Magnitude');
rmeddis@38: 
rmeddis@38: