Long time ago some clever guy emailed me this code for a transmit equaliser. In this example below I have adapted it for inverse sinc correction for simplicity. 

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

n = 31;                                  %number of taps

f = linspace(-.499,.499,n)';             %frequency range

A = 1./sinc(f);                          %required response

M = exp(-j*2*pi*f*[-(n-1)/2:(n-1)/2]); %matrix of n sinusoids 

h = inv(M'*M)*M'*A;             %scale matrix by Amp and solve equations

freqz(real(h));

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

The code looks compact. I can understand down to the line of frequency matrix (M) but no idea how is (h) coming out and it looks correct.

What is this method called? is it used in any standard filter design functions?

Any help appreciated.

The desired frequency response of h is A. That means we want DTFT(h) to be approximately equivalent to A.

The matrix M looks like a DFT matrix (or a rectangular submatrix of one). Therefore, M*h is a frequency-domain representation of h, evaluated at frequencies f. So, it looks like the idea is to find the h which gives a "good" approximation to:

M*h ≈ A.

Specifically, if we define "good" to mean optimal in a least-squares sense, then:

h = pinv(M)*A

where pinv(M) is the Moore-Penrose Pseudoinverse of M.

If M has full column rank (which it does in this case), then pinv(M) = inv(M'*M)*M'. However, we should almost never explicitly invert a matrix. Therefore, the formulation in your code snippet could suffer in terms of numerical stability and execution time. Better to just use pinv(M)*A, which I think internally uses a numerically stable Singular Value Decomposition. Or solve the system of equations directly using MATLAB's \ operator.

Thanks for the explanation.

I compared it to ifft and to my surprise I see not much difference.

So is it worth going to matrix inversion instead of ifft?

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

n = 31;                                %number of taps

f = linspace(-.499,.499,n)';           %frequency range

A = 1./sinc(f);                        %required response

M = exp(-j*2*pi*f*[-(n-1)/2:(n-1)/2]); %matrix of n sinusoids 

h = inv(M'*M)*M'*A;                    %scale matrix by Amp and solve equations

A2 = fftshift(A);

h2 = fftshift(ifft(A2));

figure;hold;

plot(real(h),'*-');

plot(real(h2),'*g--');

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

If you only care about the frequency response at n points spaced equally between DC and fsamp, then M becomes a square n-point DFT matrix.

Since DFT matrices are (scaled) unitary matrices, M'*M is then a (scaled) identity matrix, and inv(M'*M) can be replaced with a scalar 1/n.

In other words, h = inv(M'*M)*M'*A = (1/n * M') * A. Since (1/n * M') is an inverse DFT matrix (with the same 1/n scaling as MATLAB's ifft function) ... this is exactly the same as computing h = ifft(A).

If you need to do this calculation extremely fast, then this could be a useful simplification. However, I suspect there could be some value in retaining the freedom to choose the length of f independently from n. Otherwise, for small n, I could imagine large errors could occur between the DFT bins. That's just a hunch, but it could be easily investigated in MATLAB.