%
% function [Y,f,YY,y,x,xn] = dct_demo( [FS,L,AMP,FREQ[,AMP2,FREQ2[,...]]])
%
%          DCT demonstration in signal denoising
%              *** THIS IS NOT AN EXACT DENOISING ALGORITHM BUT 
%              JUST A DEMO FOR EDUCATIONAL PURPOSES ON DCT
%
%          A singal is corrupted with zero-mean random noise of uniform
%          distribution, it is then Cosine transformed and filtered.
%          Finally the signal is reconstructed using the inverse DCT
%
%          If used without parameters it implements a demo otherwise:
%          parameters are:
%              FS: the sampling frequency
%              L: the length of the signal
%              AMP: the amplitude of a simple sinusoidal
%              FREQ: the frequency of a simple sinusoidal
%              Additional parameters can be used to make the signal complex
%              The additional parameters are in pairs AMP,FREQ
%
%          output variables are:
%              Y: single-sided magnitude of the DCT 
%              f: the frequency space
%              YY: single-sided filtered magnitude of the DCT
%              y: reconstructed signal
%              x: the original sinusoid input signal
%              xn: the corrupted sinusoid input signal
%
%          Examples:
%              dft_demo;
%              [Y,f,YY,y] = dct_demo( 500, 1000, 1.2, 60, 0.8, 100);
%

%
% Copyright © 07-Nov-2012 George P. Pavlidis® (sepdek@sepdek.net, gpavlid@gmail.com)
%
function varargout = dct_demo( varargin)

	switch nargin
		case 0				% demo data
			Fs = 1000;																	% Sampling frequency
			T = 1/Fs;																		% Sample time
			L = 1000;																		% Length of signal
			t = (0:L-1)*T;															% Time vector
			x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);		% construct demo signal
		otherwise,		% user data
			if nargin<4
				error('%s ERROR: At least 4 parameters needed! See help %s',mfilename,mfilename);
			end
			Fs = varargin{1};																	% Sampling frequency
			T = 1/Fs;																					% Sample time
			L = varargin{2};																	% Length of signal
			t = (0:L-1)*T;																		% Time vector
			x = varargin{3}*sin(2*pi*varargin{4}*t);					% construct the first part of the signal
			for i=5:2:nargin,																	% for the remaining parameters
				x = x + varargin{i}*sin(2*pi*varargin{i+1}*t);	% add up more parts to the signal
			end
	end
	
	y = x + 1.5*randn(size(t));					% final noisy signal: sinusoids plus noise
	NFFT = 2^nextpow2(L);								% next power of 2 from length of y
	Y = dct(y,NFFT);										% apply the DCT
	f = Fs*linspace(0,1,NFFT);					% determine the frequency space
	
	% attempt a signal denoising by simply cutting off low amplitude
	% frequencies (below 0.25);
	YY = Y;
	dY = diff( YY);						% get derivative of YY
	stds = 5;
	%minima = find(dY<=mean(dY)+stds*std(dY) & dY>=mean(dY)-stds*std(dY));
	YY((dY<=mean(dY)+stds*std(dY) & dY>=mean(dY)-stds*std(dY)))=0;
	y_req = idct(YY,NFFT);		% the reconstructed signal after denoising by cutting off low mag freqs
	
	if nargout==0
		figure;
		% plot the noisy signal
		subplot(221); plot(Fs*t(1:100),y(1:100),'LineWidth',1);
		grid on; axis tight;
		title('Signal Corrupted with Zero-Mean Random Noise');
		xlabel('Time (milliseconds)');
		ylabel('Signal apmplitude');
		% plot single-sided amplitude spectrum
		subplot(223);  plot(f,Y,'LineWidth',1);
		grid on; axis tight;
		title('Single-Sided Amplitude Spectrum of y(t)');
		xlabel('Frequency (Hz)');
		ylabel('Magnitude |Y(f)|');
		% plot single-sided filtered amplitude spectrum
		subplot(224);  plot(f,YY,'LineWidth',1);
		grid on; axis tight;
		title('Single-Sided Amplitude Spectrum of y(t)');
		xlabel('Frequency (Hz)');
		ylabel('Magnitude |Y(f)|');
		% plot the denoised signal
		subplot(222);
		%plot(Fs*t(1:100),x(1:100),'k--','LineWidth',1);
		plot(Fs*t(1:100),y_req(1:100),'LineWidth',1);
		grid on; axis tight;
		title('Signal Denoised');
		xlabel('Time (milliseconds)');
		ylabel('Signal apmplitude');
	else
		varargout{1} = 2*abs(Y(1:NFFT/2+1));		% export the single-sided amplitude spectrum
		varargout{2} = f;												% export the frequency space
		varargout{3} = 2*abs(YY(1:NFFT/2+1));		% export the single-sided filtered amplitude spectrum
		varargout{4} = y_req;										% export the denoised signal
		varargout{5} = x;												% export the original input signal
		varargout{6} = y;												% export the corrupted input signal
	end