% Computes the power spectrum
% function [power]=  plp_ffthr(real, ll, m)
%
% real (in): Signal that is to be fft analyzed
% ll (in): number of samples of auditory spectrum               
% m (in): number of samples of auditory spectrum               
%                                                                    
% power (out): power spectrum
%

function [power]=  plp_ffthr(real, ll, m)

k=0;
n=2^m;

for ii=1:ll
	sigima(ii)=0;
	signal(ii)=real(ii);
end

for ii=ll+1:n
	sigima(ii)=0;
	signal(ii)=0;
end

n2=n/2;
n1=n-1;
j=1;
for i=1:n1,
	if i < j
		t1=signal(j);
		t2=sigima(j);
		signal(j)=signal(i);
		sigima(j)=sigima(i);
		signal(i)=t1;
		sigima(i)=t2;
	end
	k1=n2;
	while k1 < j
		j=j-k1;
		k1=k1/2;
	end
	j=j+k1;
end

for l=1:m,
	le=2^l;
	le1=le/2;
	u1=1.0;
	u2=0.0;
	w1=cos(pi/le1);
	w2=-sin(pi/le1);
	for j=1:le1,
		%for i=1:le:n
		for i=j:le:n
			id=i+le1;
			t1=signal(id)*u1-sigima(id)*u2;
			t2=sigima(id)*u1+signal(id)*u2;
			signal(id)=signal(i)-t1;
			sigima(id)=sigima(i)-t2;
			signal(i)=signal(i)+t1;
			sigima(i)=sigima(i)+t2;
		end
		u3=u1;
		if k == 0
			u1=u1*w1-u2*w2;
			u2=u2*w1+u3*w2;
		else
			u4=w1*w1+w2*w2;
			u1=(u1*w1+u2*w2)/u4;
			u2=(u2*w1-u3*w2)/u4;
		end
	end
end

if k ~= 0
	for i=1:n,
		sigima(i)=sigima(i)/n;
		signal(i)=signal(i)/n;
	end
end

for ii=1:n,
	power(ii)=signal(ii)*signal(ii)+sigima(ii)*sigima(ii);
end