% emd_2fmsin.m
%
% P. Flandrin

N = 400;
t = 1:N;

x1 = amgauss(N,N/2,fix(N/1.5)).*real(fmsin(N,.025,.2,N,3*N/4)); 
x2 = amgauss(N,N/2,fix(N/1)).*real(fmsin(N,.065,.4,N,3*N/4)*exp(i*pi/2)); 

m1 = 1.1*max(abs(x1));
m2 = 1.1*max(abs(x2));
x = x1' + x2';
m = 1.1*max(abs(x));
c = 1;
sd = .01;

figure(1)
clf

subplot(311)
plot(t,x1)
axis([1 N -m1 m1])
set(gca,'YTick',[]);set(gca,'XTick',[])
title('chirp 1','FontSize',16)
pause
subplot(312)
plot(t,x2)
axis([1 N -m2 m2])
set(gca,'YTick',[]);set(gca,'XTick',[])
title('chirp 2','FontSize',16)
pause
subplot(313)
plot(t,x,'k')
axis([1 N -m m])
set(gca,'YTick',[]);set(gca,'XTick',[])
title('chirp 1 + chirp 2','FontSize',16)

AA = axis;

pause

clf

r = x;
ner = length(x);
nzr = length(x);

imf = [];
k = 1;

while ner+nzr > 3
	
	m = r;

	nem = 1;
	nzm = length(x);
	
	mp = m;
	
	s = sd + 1;
	
	test = 0;
	
	kk = 0;

	while s > sd & test == 0,
	
		y = m;
		[indmin,indmax] = extr(y);
		
		if indmax(1) < indmin(1)
			ylmax = y(indmax(1));
			ylmin = y(1);
		else
			ylmax = y(1);
			ylmin = y(indmin(1));
		end
		
		if indmax(length(indmax)) < indmin(length(indmin))
			yrmax = y(length(y));
			yrmin = y(indmin(length(indmin)));
		else
			yrmax = y(indmax(length(indmax)));
			yrmin = y(length(y));
		end
	
		if c == 0
			envmax = interp1([t(1) t(indmax) t(length(y))],[ylmax y(indmax) yrmax],t,'linear');
			envmin = interp1([t(1) t(indmin) t(length(y))],[ylmin y(indmin) yrmin],t,'linear');
		elseif c == 1
			envmax = interp1([t(1) t(indmax) t(length(y))],[ylmax y(indmax) yrmax],t,'spline');
			envmin = interp1([t(1) t(indmin) t(length(y))],[ylmin y(indmin) yrmin],t,'spline');
		elseif c == 2
			envmax = interp1([t(1) t(indmax) t(length(y))],[ylmax y(indmax) yrmax],t,'cubic');
			envmin = interp1([t(1) t(indmin) t(length(y))],[ylmin y(indmin) yrmin],t,'cubic');
		end
	
		envmoy = (envmin+envmax)/2;
		m = y - envmoy;
		
		if k == 1 & kk <= 2
	
			subplot(2,1,1)
			plot(y,'k');hold on;
			axis(AA);
			title(['IMF ',int2str(k),';   iteration ',int2str(kk)],'FontSize',16);
			pause
			plot(t(indmax),y(indmax),'.b','MarkerSize',20)
			pause
			plot(t,envmax,'-b');
			pause
			plot(t(indmin),y(indmin),'.r','MarkerSize',20)
			pause
			plot(t,envmin,'-r');
			pause
			plot(t,envmoy,'m','LineWidth',2);
			pause
			hold off
			subplot(2,1,2);
			plot(t,m,'k');
			axis(AA);
			title('residue','FontSize',16);
			pause
		
		else

			subplot(2,1,1)
			plot(y,'k');hold on;
			axis(AA);
			title(['IMF ',int2str(k),';   iteration ',int2str(kk)],'FontSize',16);
			plot(t(indmax),y(indmax),'.b','MarkerSize',20)
			plot(t,envmax,'-b');
			plot(t(indmin),y(indmin),'.r','MarkerSize',20)
			plot(t,envmin,'-r');
			plot(t,envmoy,'m','LineWidth',2);
			hold off
			subplot(2,1,2);
			plot(t,m,'k');
			axis(AA);
			title('residue','FontSize',16);
			pause(.5)

		end
		
		s = sqrt(sum((m-mp).^2)./sum(mp.^2));
		
		[indmin,indmax] = extr(m);
		
		if length(indmin) + length(indmax) < 3
			test = 1;
		end
			
		mp = m;
		
		kk = kk + 1;
			
	end
	imf(k,:) = m;
	k = k+1;
	r = r - m;
	[indmin,indmax,indzer] = extr(r);
	ner = length(indmin) + length(indmax);
	nzr = length(indzer);
	
	pause
	
end

imf(k+1,:) = r;

f2c = [];
f2c(1,:) = imf(1,:);

c2f = [];
c2f(1,:) = imf(k+1,:);

for j = 2:k+1
	f2c(j,:) = f2c(j-1,:) + imf(j,:);
	c2f(j,:) = c2f(j-1,:) + imf(k+2-j,:);
end

emd_visuPF(x,t,imf,2);

figure(2)

pause

Nf = 256;% # of frequency bins
Nh = 81;% short-time window length
w = window(Nh,'Kaiser');

T = 1:2:N;

[s,rs] = tfrrsp(x',T,Nf,w,1);
[s,rs1] = tfrrsp(imf(1,:)',T,Nf,w,1);
[s,rs2] = tfrrsp(imf(2,:)',T,Nf,w,1);

figure(3)
clf

subplot(221)
imagesc(flipud(rs(1:128,:)))
set(gca,'YTick',[]);set(gca,'XTick',[])
xlabel('time')
ylabel('frequency')
title('signal')

pause

subplot(223)
imagesc(flipud(rs1(1:128,:)))
set(gca,'YTick',[]);set(gca,'XTick',[])
xlabel('time')
ylabel('frequency')
title('mode #1')

pause

subplot(224)
imagesc(flipud(rs2(1:128,:)))
set(gca,'YTick',[]);set(gca,'XTick',[])
xlabel('time')
ylabel('frequency')
title('mode #2')

colormap(jet)
