Eigenfilters in Matalb [closed]

I have to implement an eigenfilter for an arbitrary frequency response in MATLAB. I have this algorithm:

$N$ - order of the filter

$M = N/2$

$c(\omega) = [1,\; \cos(\omega),\; \cos(2\omega), \ldots \cos(M\omega)]'$

$d = \int D^2(\omega)$ over the desired frequency bands (Should be a scalar)

$p' = \int D(\omega)c(\omega)$ over the desired frequency bands (Should be a vector)

$Q = \int c(\omega)c(\omega)'$ over the desired frequency bands (Should be a matrix)

My objective for the process goes beyond but my concerns are limited to this part.

My first question is regarding $Q$ matrix, here I am using trapz in MATLAB but it gives me a vector instead of a matrix. I tried using loops so as to integrated is element by element (e.g. C(i*w)*c(j*w)') but then it gives me matrix of all zeros.

My second question is how to get this desired plot into a form over which I can operate in MATLAB. I can plot this using line command but I am not sure if I will be able to work on it (as I have to find $p$ vector and $d$ scalar).

I am attaching a picture for the desired frequency response. I would be really glad if anyone could help me out here. I have been banging my head with this problem for a while now.

• Are you referencing a paper or book? Is this a homework assignment? I'll bet that the paper "On The Eigenfilter Design Method and Its Applications: A Tutorial" by Tkacenko, Vaidyanathan, and Nguyen will be helpful. Is something like that what you are after? – hops Mar 27 '17 at 6:04
• It is related to that. And I have read that paper and some other paper as well. The problem is not the concept but how to implement it on matlab. As I said i have only those above questions. – Copernicus Mar 27 '17 at 6:21
• I've used this technique before. I think I could help you, if I understood the questions (and/or problems you are having) better. I've read your question as posted, but the actual question or trouble you are having is still unclear to me. You mention using trapz? Are you attempting numerical integration to obtain the matrix? When I have used this method, I generally performed the integration analytically using a piecewise description of the desired response (as you have shown), then I find the eigenvectors of the corresponding matrix equation to obtain the filters. – hops Mar 27 '17 at 7:11
• Based on jojek's edits, the first question seems to be a lack of understanding of the algorithm and of MATLAB coding. It's a little difficult to pinpoint it without more details on what you have tried exactly. The second question is a basic math question. First, you need to make sure that your units for frequency match your equation description (radians per sample in this case). Then, you can use a piecewise description of the function. I'll try to come back later when I have more time to provide a proper answer if no one else beats me to it. – hops Mar 27 '17 at 14:28

It seems to me based on your question, that you are still missing some important concepts related to the algorithm, but you seem to think that the problem is not the concept. Let me see if I can clear this up.

First of all, you seem to be referencing the book Multirate Systems and Filter Banks by Vaidyanathan (pages 53-55). This would have been helpful information that you should have included in your post (or whatever other reference you obtained this from). The exposition in the book is strictly for designing eigenfilters with a low pass response (ideally with unity gain in the passband and zero gain in the stopband). The tutorial paper On The Eigenfilter Design Method and Its Applications: A Tutorial and its references extend this method to the design of filters with an arbitrary response. You should be following those references for your problem. The key point is to use reference frequencies in the frequency response to cast the cost function in the desired quadratic form so that the Rayleigh Principle can be applied to find a minimum solution (namely, the eigenvector corresponding to minimum eigenvalue).

I'd be happy to help walk you through my understanding of the filter design procedure for your particular case, but I have limited free time. For this reason, I won't do so unless I hear that you are interested.

Okay, since this is a Q&A website, and I want to keep this short, I will be glossing over a lot of details that you should definitely read on your own. The tutorial paper referenced above indicates that we can obtain an arbitrary response for the eigenfilter by minimizing the cost function $$\xi = \mathbf{h}^\dagger \mathbf{P} \mathbf{h}$$ where I have adhered to the notation in the paper so $\dagger$ indicates the conjugate transpose (a.k.a. hermitian). Ignoring the constants and the weighting function from the paper, the matrix $\mathbf{P}$ can be computed as $$\mathbf{P} = \int_{\mathcal{R}} \left[ \frac{D(\omega)}{D(\omega_0)} \mathbf{e}(e^{j\omega_0}) - \mathbf{e}(e^{j\omega}) \right]^* \left[ \frac{D(\omega)}{D(\omega_0)} \mathbf{e}(e^{j\omega_0}) - \mathbf{e}(e^{j\omega}) \right]^T d\omega$$ where $\mathcal{R}$ is the region over which $D(\omega)$ is defined, $D(\omega)$ is your desired frequency response, $\omega_0$ is a fixed reference frequency, and $\mathbf{e}(z) = \left[\begin{array}{cccc} 1 & z^{-1} & \cdots & z^{-N} \end{array} \right]^T$. When I have more time, I will provide more details on how to go about computing this (and converting it from being complex-valued to real-valued and linear phase). Although, this is computable using trapz as you suggested in your original question, and so it may be an acceptable solution already.

I believe that your specific case for $D(\omega)$ can be defined as follows $$D(\omega) = \left\{\begin{array}{cc} 1 & |\omega| \leq \frac{\pi}{3}\\ 0 & \frac{5\pi}{12} \leq |\omega| \leq \frac{7\pi}{12} \mbox{ or } \frac{9\pi}{12} \leq |\omega| \leq \pi\\ \frac{12}{\pi}\left(|\omega| - \frac{15 \pi}{24}\right) & \frac{15\pi}{24} \leq |\omega| \leq \frac{17 \pi}{24} \\ \mbox{undefined} & \mbox{otherwise} \end{array}\right.$$ where $\omega$ can vary between $-\pi$ and $\pi$. Let me know if that is correct. This gives you a total of 7 nonzero regions over which you need to integrate to obtain the elements of the matrix $\mathbf{P}$.

A convenient reference frequency for you will be $\omega_0 = 0$ or DC. This leads to a special form where $\mathbf{e}(e^{j \omega_0}) = \mathbf{1}$ and $\mathbf{1}$ is a column vector of all ones with the same length. Also, $D(\omega_0) = 1$ so that simplifies things a bit. For now, if you choose a spectrum with conjugate symmetry, you should obtain a real response. It would be best to re-derive the formula in terms of $\mathbf{c}(\omega)$ (a vector of cosines analogous to the paper's definition for $\mathbf{e}(\omega)$ that exploits the symmetry in the coefficients). You might attempt to use this information together with the equation for the matrix $\mathbf{P}$ to find your filter coefficients. If you do and post the MATLAB that you generate, then I can do more to assist you with this.

• Sorry for the late reply, I was out of town. I would be glad if you could help. I have read the paper and it talks about the reference frequency but the algorithm I am using is almost like that but not exactly. It doesnt require a reference frequency. It could be the lack of understanding the algorithm and matlab that I am not able to do it. I would be glad if you help me out here. – Copernicus Mar 28 '17 at 18:08
• I'll attempt to provide a framework for your particular case when I have some free time. Hopefully this week sometime. – hops Mar 29 '17 at 4:32
• Thank you @hops, I am half way through it now and will be done today probably. I really appreciate your efforts here. – Copernicus Apr 6 '17 at 19:52
• Hello @hops, I was able to find the matrix and everything but the response is not even close to what I want. I believe the problem comes when I find the desired frequency response of the linear part of the plot (second pass-band). Following is the code snippet where I find the the desirable response and a corresponding vector after using 'trapz': p3 = zeros(1,(N/2)+1); x = (1.9635:6.2782e-04:2.2253); y = (3.8462.*x)-7.5577; d3 = trapz(y.^2); for i = 0:N/2 p3(i+1) = trapz(w3,y.*cos(i*w3)); end – Copernicus Apr 7 '17 at 19:24
• I used your equation for the band and it gives the same result! – Copernicus Apr 7 '17 at 19:27