Hi, 

I'm currently trying to implement the FFT-based equation method for filter design outlined by JOS, details of which can be found here. Specifically, I'm having trouble with how to implement the equations using Toeplitz matrices and from which the b, a coefficiants are derived. 


For the above expression it isn't clear (to me) exactly how the indexing would work for R_yu.

I'm not also not entirely sure how the below expression operations. Am I right in thinking that it's some kind of nested matrix multiplication e.g R_uu matrix multiplied with R_yu, the result of which is then added to the matrix multiplied result of R_uy and R_yy? 

Regards,

Robert

Apologies, I had tried to copy and paste an image in but that doesn't seem to have worked. 

If you follow the link provided in the original post, it is the indexing involved in the generation of the Toeplitz matrices that I'm struggling to understand. And then the subsequent matrix concatenation operation which results in the final theta parameters, containing the ceofficients. 

Obviously Jos is the best to reply. Admittedly the math in his document is well above my capacity. However, I have come across some mystery matrix method of FIR filter design instead of iDFT that may be related to your post. In this method you generate a square matrix of sinusoids across both frequency grid and time domain then scale each and solve the equations. I have used it for equaliser FIR with some success over iDFT based iterative methods.

This is example of what I concluded but can't pretend I know the details fully.

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

n = 81;                                 %number of taps

fx = linspace(-.4999,.4999,n)';         %frequency grid

Amp = sinc(fx);                         %example required response

F = exp(-j*2*pi*fx*[-(n-1)/2:(n-1)/2]); %square grid of sinusoids

a = inv(F'* F)*F'*Amp;                  %scale and solve matrix

%%%%%%%%%% check with iDFT method %%%%%%%%%%

aa = ifft(Amp);

t1 = abs(ifftshift(fft(aa)));

t2  = abs(fft(a));

plot(t1,'*-');hold;

plot(t2,'g*--');

Hmm... by linked page you're trying to implement Matlab's invfreqz() or stmcb() ? 

Here are couple link to implementations in Python language as an example for you:

https://github.com/yurochka/dsptools/blob/master/d...

https://github.com/nbara/python-meegkit/blob/maste...