view userProgramsTim/mellin_trafo.m @ 38:c2204b18f4a2 tip

End nov big change

function out = mellin_trafo(inx,iny)

%This function computes the Mellin Transformation of a one-dimensional
%Signal in analytical terms given as
%
% S(p)=Integral from 0 to infty of s(t) t^(p-1)dt
% Which equals 
% S(c) = Integral from -Infty to Infty s(t)*exp(j*c*log(t))d(log(t))
%
% taken from Irino and Patterson, Speech Communication 2002 Eq. (1) and (3)
% Tim Juergens, September 2011
%

%Resample the Signal onto a log(t) axis
minimalx=min(inx);
maximalx=max(inx);
logarithmicx= exp([log(minimalx):(log(maximalx)-log(minimalx))/length(inx):log(maximalx)]);
logarithmicy= interp1(inx,iny,logarithmicx,'linear','extrap');
%figure, semilogx(logarithmicx,logarithmicy);

%Absolute of the inverse Fourier-Transform of the resampled signal
out = abs(ifft(logarithmicy));
%figure, plot(out(1:40));


%%just to show that the mellin transform results in invariable patterns if
%%the formants are a constant ratio
% frequencies_short = [1:8:8000];
% frequencies_long = [1:10:10000];
% Intervals_long = 1./frequencies_long;
% Intervals_short = 1./frequencies_short;
% contoursin = sin(2*pi*0.00015.*frequencies_long).^2
% figure, plot(frequencies_long,contoursin), hold on, plot(frequencies_short,contoursin,'r')
% xlabel('Frequency (Hz)')
% ylabel('Rate (arb. units)')
% figure, semilogx(Intervals_long,contoursin), hold on, plot(Intervals_short,contoursin,'r')
% xlabel('Interval (s)')
% ylabel('Rate (arb. units)')
% m1 = mellin_trafo(Intervals_long,contoursin);
% m2 = mellin_trafo(Intervals_short,contoursin);
% figure, plot(m1(1:40)), hold on, plot(m2(1:40),'r');
% xlabel('Mellin variable c');
% ylabel('Magnitude');