LAB 8, 03/26/14
________________
		

1.	Signal approximation plays an important role in audio and image 
	compression. 

	Consider the problem of compressing a vector s by a factor of 4, 
	e.g., from length N = 256 to N/4 = 64.  For example,
		
		N = 256 ;		
		t = 0:N-1 ;

		s = [((0:127)/128).^10   -((127:-1:0)/128).^10].' ;
		
		plot(t,s), axis([0 255  -1.2 1.2])

	Suppose that we keep every fourth sample in s, then interpolate 
	s(1:4:256) in the time domain using

		s1: a cubic spline

		s2: the continuous time signal of lowest bandwidth that
		    interpolates s(1:4:256) in a periodic fashion
	
	These signals are obtained as follows:

		q = 0:4:255 ;

		s1 = interp1 ( q, s(q+1), t, 'spline' ) ;
                                   % zero indices not allowed!

		F = fft ( s(q+1) );        % again!
		s2 = 4*ifft([F(1:32);F(33)/2;zeros(191,1);F(33)/2;F(34:64)]);

		subplot(2,1,1), plot(t,s,'r',t,s1),
		title(['Cubic Spline Interpolation']),
		axis([0 N-1 -1.2 1.2])

		subplot(2,1,2), plot(t,s,'r',t,s2),
		title(['Bandlimited Shannon-Nyquist Interpolation']),
		axis([0 N-1 -1.2 1.2])

	Although s1 is a better interpolant than s2, the agreement is not
	particularly good near the spike (in the graph of s).


2.	Alternatively, we can project s onto the columns of an orthonormal
	NxN matrix, keep the N/4 projection coefficients with the largest 
	absolute value, and approximate s using the corresponding matrix 
	columns.  Two such approximations are considered here:

		s3: based on the discrete Fourier transform

		s4: based on the Haar wavelet developed in Lab 7

	The results are as follows:

		S = fft(s) ;
		[temp, J] = sort(abs(S),'descend') ;
		J = J(1:N/4) ;
		S = S.*(ismember((1:N)',J)) ;
		s3 = real(ifft(S));

		V = nhaar222(N) ;
		c = V'*s ;
		cabs = abs(c);
		[temp, J] = sort(cabs,'descend') ;
		J = J(1:N/4) ;
		c = c.*(ismember((1:N)',J)) ;
		s4 = V*c ;

		subplot(2,1,1), plot(t,s,'r',t,s3),
		title(['Approximation By DFT Compression']),
		axis([0 N-1 -1.2 1.2])

		subplot(2,1,2), plot(t,s,'r',t,s4),
		title(['Approximation By Haar Wavelet Compression']),
		axis([0 N-1 -1.2 1.2])		

	For this particular signal, the Haar Wavelet method gives by far the
	most satisfactory result; i.e., s4 is a better approximation to s 
	than any of s1, s2 and s3.

	
3.	We now examine the syntax of user-defined functions, such as 
	HAAR222 and NHAAR222 created in the last lab (also used above in 
	creating the matrix V).  We display NHAAR222, namely the contents of 
	nhaar222.m, below:

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	function V = nhaar222(N)

	% NHAAR222 Normalized Haar wavelet matrix.
	% 	V = NHAAR222(N), where N is a positive integer, is an 
	%	orthonormal square matrix of size MxM, where M is the largest 
	%	integer power of 2 less than or equal to N. The first column 
	%	of the matrix is constant, while the remaining columns are
	%	Haar wavelets at log2(M)-1 different scales.


		K = fix(log2(N)) ;
		a = ones(N/2,1) ;
		v0 = [a ; a];
		v1 = [a ; -a];
		V = [v0 v1];
		for i=2:K
			v = v1(1:2^(i-1):N);
			v = [v ; zeros(N-length(v),1)];
			V = [V v];
			for j = 1:2^(i-1)-1		
				v = circshift(v,2^(K-i+1));
				V = [V v];
			end
		end

		V = V*diag(1./sqrt(diag(V'*V)));

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

4.	- The first line is the function definition line, always beginning
	   with "function" in lower-case letters.  The statement

			V = nhaar222(N)

	  gives the name of the function (NHAAR222), the number of input
	  arguments (1) and the number of output arguments (1).  The variables
	  N and V are local to the function.  Although the names N and V are
	  arbitrary, their use must be consistent throughout the body of the
	  function.

	- The first comment line (line 'H1') is a brief description of the
	  function, that comes up when using the LOOKFOR search command:

		lookfor nhaar222

	  The H1 line together with the comment lines below it provide the
	  on-line documentation of the function, displayed using HELP:

		help nhaar222


5.	- The body of the function is a set of commands executed in succession
	  (i.e., as a script).  Briefly, V is created as follows:

		- N is truncated to 2^K, using the function FIX
		  (related functions: ROUND, CEIL, FLOOR);

		- the first two columns (Scale 0 and Scale 1) are formed;

		- for each scale i=2:K, the Scale 1 column is compressed by
		  sampling every (2^(i-1))th entry, then zero-padded to
		  produce the first column at Scale i;

		- the first column at Scale i is circularly shifted 2^(i-1)-1
		  times to produce the remaining columns at that scale;
		
		- each column is normalized by dividing it by its norm.


6.	User-defined functions:

		- can have any number (including zero) of input and
		  output arguments

		- can execute different tasks depending on the number
		  of input and output arguments (this is achieved through 
		  conditional statements using NARGIN and NARGOUT)

	In the earlier example (2. above), where the signal s was compressed
	by a factor of 4 using the N/4 (=64) absolutely largest coefficients,
	we used the code

		c = V'*s ;
		cabs = abs(c);
		[temp, J] = sort(cabs,'descend') ;
		J = J(1:N/4) ;
		c = c.*(ismember((1:N)',J)) ;
		s4 = V*c ;

	Note that in the syntax used above, the function SORT has 

	- two input arguments: the vector cabs and the string 'descend'; and

	- two output arguments: the sorted version of cabs (denoted by temp)
	  and the index vector J such that temp = cabs(J).

	The first N/4 indices in J are kept, and each of 1:N is tested for
	inclusion in the vector J using the ISMEMBER function:

		ismember((1:N)',J)

	returns a logical 0-1 array of the same size (N). Multiplying each
	entry of c by the corresponding entry of that array results in nulling
	out the insignificant coefficients of c.


7.	Our final example is related to the script in 6. above.  We write
	a function COMPRESS1 that finds the M absolutely largest entries
	of a real or complex vector X, nulls out the remaining entries, 
	and also computes the energy (square norm) ratio between the output 
	and input vectors.


	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	function [y,J,r] = compress1(x,M)

	% COMPRESS1 Fixed-length compression by eliminating smallest entries.
	% 	[Y,J,R] = COMPRESS1(X,M) determines the indices J containing 
	%	the M absolutely largest entries of a vector X, nulls out the 
	%	remaining entries and also computes the square norm ratio 
	%	between the output vector Y and the input vector X.

		M = min(M,length(x)) ;
		xabs = abs(x) ;
		[temp, J] = sort(xabs,'descend') ;
		J = J(1:M) ;
		I = reshape(1:length(x),size(x)) ;
		y = x.*(ismember(I,J)) ;
		r = (norm(y)/norm(x))^2 ;

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%