2
$\begingroup$

I'm studying FIR filter design and it's time for the frequency sampling method, my teacher said that to use this method you need to follow the following steps:
  • Sample the periodic frequency response of the ideal filter you want to approximate;
  • Calculate the IDFT of the samples;
  • Use the results as filter coefficients.

I tried to design a lowpass filter on MATLAB and it seems to be working but I can't figure out why this method works this way, and what it is practically doing. Why do I need to periodize and then sample the frequency response over one period? Since the sinc() function is not limited in time should I get some aliasing in the time domain while sampling the frequency response? What I'm trying to say is that I understood how to write the code but not how this method works behind the scenes, I'm trying to extract a sort of step-by-step guide, about how this type of FIR design works according to signal processing theory, like the following (it may be wrong, I don't know):

  1. Sample the frequency response of the ideal filter (just one rect centred in f=0Hz in case of lowpass);
  2. the sampling will give you the periodization of the sinc() in the time domain;
  3. sampling the sinc() you'll have the periodization of the frequency response;
  4. and so on ...

because I'm sure it will help me to better understand this method and to write it down in my notebook, I'm interested just in one theoretical explanation. I hope you can help me.

$\endgroup$
1
  • $\begingroup$ "Sample the periodic frequency response of the ideal filter you want to approximate...". Are you starting with a sampled-time signal, and that's why it's periodic? I know this is a nit (and I'm going to assume that you are, in fact, working in the sampled-time domain all along), but please edit your question to say whether you're starting out with a sampled-time but infinite-length prototype, or a continuous-time prototype. $\endgroup$
    – TimWescott
    Commented Apr 26 at 16:19

5 Answers 5

3
$\begingroup$

Explanation:

The impulse response and frequency response of a filter are Fourier Transform pairs.

The time domain response of a filter is its impulse response: you put an impulse at the input of the filter and the response from that impulse will appear at the filter's output as time moves forward.

The Fourier Transform of an impulse is all frequencies, so by applying an impulse at the input in the time domain, we are equivalently applying all frequencies in the frequency domain.

The frequency domain output of the filter with all frequencies applied at the input is then the Frequency Response of the filter.

For Digital Filters, FIR filters specifically, the coefficients of the filter are samples of the impulse response. Sure enough, if you apply the unit sample (discrete time impulse) at the input of an FIR filter, it will respond with the coefficients at the output as that unit sample moves through the filter.

coefficients of FIR Filter

Thus by having samples of the desired impulse response, we can use that as the filter coefficients for an FIR filter.

Frequency Response from Coefficients

The FFT is an algorithm that provides the Discrete Fourier Transform (DFT): this would be samples of the Frequency Response. (A detail here is the Frequency Response for a discrete-time filter such as an FIR is continuous in the frequency domain, and is specifically the Discrete-Time Fourier Transform or DTFT. The DFT is samples of the DTFT and therefore samples of the frequency response).

Thus the inverse FFT of a target frequency response provides us with samples of the impulse response (almost). There is a distortion involved in doing this called Time Domain Aliasing, so the samples will not exactly be the impulse response, but an "aliased" impulse response. If the time duration given by the number of samples used is shorter than the actual impulse response needed for the target frequency response, then those tails that go beyond that time duration alias back into the time duration given.

Increasing the total time duration (by oversampling in the frequency domain: increasing the number of bins in the FFT) reduces this aliasing in the time domain, and thus results in a more accurate estimate of the actual samples for the impulse response that would result in the frequency response as given by the FFT. (See Robert's answer discussing how to optimize the filter design for this- which I refer to as the Windowing Method which this becomes identical to when we sufficiently eliminate the time aliasing distortion, which is a different and typically better frequency design approach).

The pure Frequency Sampling approach as the OP has described is NOT a very effective method for designing filters, and hopefully this is made clear when it is introduced. It is very simple and useful when we need an exact solution at specific frequencies only (such as in OFDM modulation), but for a continuous frequency response it will have more error than other approaches given the same complexity (same number of taps). Specifically the time domain aliasing mentioned above results in greater pass band error and less stop band rejection than we could otherwise get with approach such as Windowing that RBJ mentions in his answer. I compare these approaches directly as well as others in DSP.SE #31905. The plot below shows specifically the time domain aliasing that occurs when using these two approaches: the orange curve is a pure Sinc function, as the inverse Fourier Transform of the rectangular frequency response as a target frequency response. The blue curve shows what we would get if we used an inverse FFT of samples of that rectangular response given the number of samples used (time domain aliasing!). The pure Sinc as shown is only a portion of the ideal Sinc that extends to +/-infinity, so distortion due to truncation will result (and hence the use of windowing to reduce that). Still, the distortion due to truncation of the pure Sinc will be less than the distortion from the time domain aliasing, as I demonstrated in more detail at the linked question mentioned above.

Frequency Sampling

The approach of using an increased time duration (increasing the total number of samples in the DFT which results in a longer time duration inverse FFT given the same sampling rate), meaning significantly oversample the desired frequency response to result in a sufficiently accurate time domain impulse response (minimized time domain aliasing), and then windowing this to a smaller filter size, is in contrast a very effective and accurate design approach, and can certainly be referred to as "Frequency Sampling" since the desired frequency response is sampled. However, this is more consistent with the "Windowing Method" where the desired accurate impulse response (however that is obtained) is sampled and then windowed and I think that distinction as well as the poor results with pure "Frequency Sampling", is important to understand.

$\endgroup$
2
  • $\begingroup$ Thanks again, I think I almost got it, does the "explanation" I posted as an answer makes any sense? $\endgroup$
    – minghierid
    Commented Apr 27 at 9:05
  • 1
    $\begingroup$ @minghierid yes it did, nice drawings too! $\endgroup$ Commented Apr 28 at 20:18
4
$\begingroup$

This is how I would define the frequency-sampling method.

Assuming linear phase (to start with).

  1. Draw out your desired frequency response magnitude, with as much detail as you need to to define the features you want. This is the "target frequency response".

  2. Put that in as the real part of a huge FFT. Like 1024K in size. This means that Nyquist/512K is your frequency resolution (let's hope it's a lot more than you will need). Reflect the "negative frequencies" to the upper half of the FFT.

  3. Perform the inverse FFT. Now you will have, centered on the 0th bin (the latter half of the result are your negative-time samples), an impulse response, but it's 1024K samples long. Too long to be practical.

  4. Apply a Kaiser window to that too long impulse response to make it shorter. Center this window at the 0th bin. Zero all samples outside the window width. You still have 1024K samples, it's just that most of them are zero. And it's still a symmetrical impulse response.

  5. FFT that result to see what frequency response you're gonna get.

  6. If you don't like that frequency response, let's do this again, but make the Kaiser window wider. Or adjust the $\beta$ parameter in the Kaiser window. Se what you're gonna get. Increase the length of the window (and impulse response) until the resulting frequency response is within acceptable limits of deviation from your original target.

That is, according to me, how the frequency-sampling FIR method works.

$\endgroup$
1
  • $\begingroup$ Minor detail but I actually refer to this as the windowing method of filter design, as that method is specifically windowing the time domain impulse response, which you have done. Here an oversampled inverse FFT was used to estimate the impulse response while minimizing the time domain aliasing we would otherwise get with the Frequency Sampling method directly. What the OP describes is what I consider the "Frequency Sampling Method" it its purest form. (Your method is much better, I just make this distinction from the Windowing Method). $\endgroup$ Commented Apr 26 at 3:18
4
$\begingroup$

This answer shows one way of understanding why and how frequency sampling works. Let's assume you're given a desired complex frequency response $H(\omega)$ on the interval $\omega\in[0,\pi]$, where $\omega=\pi$ corresponds to the Nyquist frequency. First, we need to extend the frequency response to the interval $[0,2\pi]$, in order to have one full period of the periodic spectrum. Since we want to design a filter with real-valued filter coefficients, we have the following symmetry:

$$H(\omega)=H^*(2\pi-\omega)\tag{1}$$

Assuming that $H(\omega)$ is specified for all $\omega$ and using $(1)$, we could find the corresponding impulse response by computing the inverse discrete-time Fourier transform (IDTFT):

$$h[n]=\frac{1}{2\pi}\int_0^{2\pi}H(\omega)e^{jn\omega}d\omega\tag{2}$$

In general, $h[n]$ will be infinitely long and non-causal.

Now let's approximate the integral $(1)$ by a Riemann sum. We use a uniform frequency grid

$$\omega_k=\frac{2\pi k}{N},\qquad k=0,1,\ldots,N-1\tag{3}$$

with frequency distance $\Delta\omega=2\pi/N$. The corresponding Riemann sum then becomes

\begin{align*} h[n]&\approx\frac{1}{2\pi}\sum_{k=0}^{N-1}H(\omega_k)e^{jn\omega_k}\Delta\omega\\ &=\frac{1}{2\pi}\sum_{k=0}^{N-1}H(\omega_k)e^{j2\pi nk/N}\frac{2\pi}{N}\\ &=\frac{1}{N}\sum_{k=0}^{N-1}H[k]e^{j2\pi nk/N},\qquad n=0,1,\ldots,N-1\tag{4} \end{align*}

where I've used the notation

$$H[k]\triangleq H(\omega_k)$$

Note that $(4)$ is just the inverse discrete Fourier transform (IDFT) of the sampled frequency response $H[k]$.

Approximating $(2)$ by $(4)$ means that instead of a possibly infinitely long sequence $h[n]$, we obtain a sequence of length $N$. If we tried to evaluate $(4)$ for values $n$ outside the interval $[0,N-1]$, we would see that the sequence obtained by $(4)$ is periodic with period $N$. This periodization is a consequence of sampling the frequency response, and it is referred to as time-domain aliasing.

The approximation of $(2)$ by $(4)$ is the frequency sampling method for the design of FIR filters. At the frequencies $\omega_k$, the resulting frequency response equals the desired values $H[k]$. In order to make the approximation error small over the continuous frequency range, we would usually choose a large value $N$, i.e., many frequency points $\omega_k$. This will result in large filter lengths which may not be practical. That's why for most practical purposes we would use a window function to truncate the filter's impulse response of length $N$ to some lower value, which will result in a lower computational complexity, and - sometimes even more important - in a lower overall delay.

Frequency sampling could be used if explicitly solving the integral $(2)$ is not an option. This is the case if

  1. no analytic expression is given for $H(\omega)$, for instance because the frequency response is obtained by measurements at discrete frequencies,
  2. there is no closed-form solution of the integral for the given frequency response.

In these cases, frequency sampling can be a practical solution. In my experience, frequency sampling is especially useful if the given frequency response is relatively smooth and/or if the desired phase response is non-linear. If combined with windowing, it will be very close in performance to a least squares solution at a much lower computational cost.

As a final note, I don't think that frequency sampling is a useful method for designing frequency-selective filters with piecewise constant desired responses, because in that case the integral $(2)$ can be solve analytically and the well-known window method could be used, if design complexity is an issue. Otherwise, for linear-phase FIR filters, there are several optimal design methods resulting in better filters than the ones designed by frequency sampling or windowing.

In sum, the frequency sampling method for designing FIR filters should be considered if one or more of the following points are satisfied:

  • design complexity plays a role, such as in certain real-time applications
  • the desired frequency response is not one of the standard responses, and it is relatively smooth
  • the desired phase response is non-linear
  • the desired response is given by an array of (measured) values instead of by an analytical expression
$\endgroup$
8
  • $\begingroup$ Very nice and concise. Do you recall the case we discussed that I came across where if the frequency response was increasing as we approach DC (I was modelling phase noise curves) that even though it was relatively smooth oversampling the response in frequency did not help. I'm still interested in why beyond an intuition I have ---if you have further insights maybe I'll pull that together again and post it as a question. $\endgroup$ Commented Apr 26 at 21:16
  • $\begingroup$ @DanBoschen: Thanks Dan! Could you share the link to that question you refer to, that would help refresh my memory ... $\endgroup$
    – Matt L.
    Commented Apr 26 at 21:47
  • $\begingroup$ I'll dig - not sure if it was a comment in a question (which I'll have difficulty finding) or an email- you shared some initial insight so would be good to find. I need to refresh my memory as well. I'll probably recreate the case and post it. $\endgroup$ Commented Apr 26 at 23:16
  • $\begingroup$ ! Amazingly I found it: dsp.stackexchange.com/a/88072/21048. I did mention in the comments my suspicion on smoothness being involved and not really sure now what exactly I was getting at (perhaps minimum time domain response duration): as in that case the interesting point was oversampling and truncation had no effect. I have some other insights but need to solidify them. $\endgroup$ Commented Apr 26 at 23:35
  • $\begingroup$ Is the $(2)$ the IDTFT or the IDFT ? $\endgroup$
    – minghierid
    Commented Apr 26 at 23:58
2
$\begingroup$

Interesting to see the links and history from Dan and Robert and reference to Matt's thesis and derivation for FDLS, which I hadn't viewed previously. This is similar to approaches we've used since the 90's in various FIR filter optimizations. Though my own derivation proceeds slightly differently, I think it ends up in essentially the same place. The frequency sampling approach is then a degenerative subset.

Starting with a desired, complex response $G(e^{j\theta})$, we want to find $H(z)=h_0+h_1z^{-1}+\cdots+h_pz^{-p}$ as an order $p$ FIR to minimize the $L^2$ norm as \begin{eqnarray} {\bf h_{LS}} &=& arg\min_{\bf h}\int_{-\pi}^\pi \left|G(e^{j\theta})-H(e^{j\theta})\right|^2d\theta\\ &=& arg\min_{\bf h}\int_{-\pi}^\pi \left|G(e^{j\theta})-\psi^H(e^{j\theta}){\bf h}\right|^2d\theta \end{eqnarray} where $$ \psi^H(e^{j\theta}) = \left[1 \;\;e^{-j\theta}\;\;e^{-j2\theta}\;\cdots\;e^{-jp\theta} \right] $$ $$ {\bf h} = \left[ h_0 \;\; h_1 \;\; h_2 \; \cdots \; h_p\right]^T $$ We often want to have a spectral weight factor, either to step on the scale a bit at certain frequencies (e.g. to overweight a stopband rejection more heavily than in-band ripple) or to give proper spectral weight to compute the resulting EVM metric for an expected pulse shape and resulting PSD (e.g. UMTS RRC shape and EVM). This adds a design-specified spectral weight of $|W(e^{j\theta})|^2$ to the $L^2$ error minimization as $$ {\bf h_{LS}} = arg\min_{\bf h}\int_{-\pi}^\pi \left|W(e^{j\theta})\left[G(e^{j\theta})-\psi^H(e^{j\theta}){\bf h}\right]\right|^2d\theta $$ Up to this point, $W(e^{j\theta})$ is a user-defined spectral weight option (can use unity if no weight variation is desired), $G(e^{j\theta})$ is the complex, desired frequency response, and $\psi(e^{j\theta})$ is the Vandermonde column used to compute the DTFT point of the FIR coefficients in $\bf h$ at $\theta$.

Continuing, write the integral as a Riemann sum, ignoring scaling constants that don't impact the least squares solution: \begin{eqnarray} {\bf h_{LS}} &=& arg\min_{\bf h}\int_{-\pi}^\pi \left|W(e^{j\theta})\left[G(e^{j\theta})-\psi^H(e^{j\theta}){\bf h}\right]\right|^2d\theta\\ &=& arg\min_{\bf h} \lim_{N\rightarrow\infty}\sum_{k=-N/2}^{N/2-1}\left|W(e^{j\theta_k})\left[G(e^{j\theta_k})-\psi^H(e^{j\theta_k}){\bf h}\right]\right|^2\;\;\;;\;\;\theta_k = 2\pi k/N \\ &\approx& arg\min_{\bf h} \sum_{k=-N/2}^{N/2-1}\left|W(e^{j\theta_k})\left[G(e^{j\theta_k})-\psi^H(e^{j\theta_k}){\bf h}\right]\right|^2\;\;\;\;\mbox{for large } N \\ &=&\left||{\bf W}\left(\gamma - \Psi{\bf h} \right) \right||_2^2 \end{eqnarray} where $$ {\bf W} = diag\left[\left\{W(e^{j\theta_k})\right\}_{k=-N/2}^{N/2-1}\right] $$ $$ \gamma = \left[ \begin{array}{c} G(e^{-j\pi}) \\ G(e^{-j2\pi(N/2-1)/N}) \\ \vdots \\ G(e^{-j2\pi/N}) \\ G(e^{j0}) \\ G(e^{j2\pi/N}) \\ \vdots \\ G(e^{j2\pi(N/2-1)/N}) \end{array} \right] $$ $$ \Psi = \left[\begin{array}{c} \psi^H(e^{-j\pi}) \\ \psi^H(e^{-j2\pi(N/2-1)/N}) \\ \vdots \\ \psi^H(e^{-j2\pi/N}) \\ \psi^H(e^{j0}) \\ \psi^H(e^{j2\pi/N}) \\ \vdots \\ \psi^H(e^{j2\pi(N/2-1)/N}) \end{array} \right] $$ Here, $\bf W$ is the diagonal, spectral weight matrix, $\gamma$ is the vector of complex, desired frequency response values, and $\Psi$ is the (Hermetian transpose) Vandermonde matrix for computing the DTFT of $\bf h$ at N discrete, oversampled points.

The weighted LS solution is then $$ {\bf h}_{LS} = ({\bf W}\Psi)^+{\bf W}\gamma $$ where $({\bf W}\Psi)^+$ is the pseudo-inverse of ${\bf W}\Psi$.

In practice, I find that $N$ being 1 to 2 orders of magnitude larger than $p$ works well to avoid under-specifying the spectral response.

I believe the "frequency-sampling" approach is the special case of ${\bf W} = {\bf I}$ and $N=p+1$ so that $\Psi$ is the DFT matrix - within a row rotation of the $[0, 2pi)$ variant of the "DFT matrix" as it is more commonly written.

$\endgroup$
3
  • $\begingroup$ Nice vml, it is interesting how the DFT is just another least squares solution which if I follow correctly is the conclusion you reached here. Please see this related question I had posted: dsp.stackexchange.com/q/87326/21048 which is about that along with an answer I added for further detail showing how the DFT solution is inferior (consistent with "Frequency Sampling" being inferior to "Least Squares" for filter design). You may have some additional insights to that question as posted, and am still seeking an answer with some practical insights. Thanks! $\endgroup$ Commented Apr 26 at 10:20
  • $\begingroup$ Thanks, Dan, I'll have a look this weekend. As noted, the degenerate version of the above equations is a row rotation of the DFT matrix (previously mentioned incorrectly as "column rotation" in my quick write-up last night). In that degenerate case where the overfit LS setup becomes square with an exact solution, we have the simple inverse DFT answer that allows the DTFT between DFT points to move away from desired without explicit metric cost. $\endgroup$
    – vml
    Commented Apr 26 at 13:59
  • $\begingroup$ Ah very cool, I think I see what you're getting at in that it least squares without being overdetermined-- so no effective averaging of the errors (perhaps)....you have good insight into this. I'm looking forward to your perspective on that post! Thank you $\endgroup$ Commented Apr 26 at 17:33
2
$\begingroup$

Thanks to all for your help! Here's what I understand, please correct me if I'm wrong
When performing a time sampling using $T$ as the sampling step, it produces a frequency periodicization with $1/T=fc$ as the frequency step, which is known as the DTFT. Now, if we sample the periodic frequency spectrum with $\Delta f$ as the frequency step, we get a periodicization of the time samples with step $1/\Delta f$. If the samples in the time domain are in a finite number (for example, N samples from 0 to N-1), then using $\Delta f = 1/NT$ as the frequency sampling step avoids time domain aliasing, which is known as the DFT. enter image description here

Let's suppose I want to design a low-pass filter with the frequency sampling method.

enter image description here

If I only sample the frequency response, the time domain signal will be periodic with some aliasing. (does the DTFT works also this way? so sampling in frequency equals periodization in time?) This happens because the sinc() function is not time-limited. In other words, the signal will be the sum of infinite continuous sinc() functions in the time domain, resulting in another infinite duration signal.

enter image description here

I am looking to calculate the coefficients of the filter. To do this, I need to sample the impulse response (the sum of all the infinite sinc) which is now periodic (each period has a finite length), however, sampling the impulse response involves periodizing the frequency spectrum, which makes the time and frequency domain periodic. Therefore, I can consider only one period, which means N samples.

enter image description here

Since both the time and frequency domains are periodic, I can consider just one frequency response period and obtain the time domain period using IDFT on the samples (sampled impulse response). Now I can also apply any windowing function to the filter coefficients to partially get rid of that aliasing and improve the filter response Does it makes any sense?

$\endgroup$
6
  • 2
    $\begingroup$ It's true that sampling in one domain causes periodization in the other domain. In general, periodization causes aliasing, which can occur in the frequency domain as well as in the time domain. What I don't understand is why you start out in the continuous time domain. If we're talking about FIR filter design, we are already in the discrete time domain, so we can use the DTFT. In the case of frequency sampling, we also apply sampling in the frequency domain which causes time domain aliasing, but which allows us to use the DFT. The latter is efficiently implemented by the FFT. $\endgroup$
    – Matt L.
    Commented Apr 27 at 8:55
  • 1
    $\begingroup$ yes this all looks good and I love your hand written drawings. I am not as concerned as Matt is with your use of continuous time models if they help your intuition (it helps me and builds an understanding of sampling in general!), as they are valid mathematical constructions. Just to clarify your last comment on windowing, to be clear, windowing will reduce the sidelobes in the Freq Domain result, so yes less aliasing due to that. However notice the techniques described here are to increase the samples in frequency first (extend the time domain) and then truncate/window. $\endgroup$ Commented Apr 27 at 11:45
  • 1
    $\begingroup$ So what we are doing with that is getting a better approximation of the optimum impulse function we would have without any aliasing, before any windowing is applied! The Sinc vs the Dirichlet Kernel (aliased Sinc) is a great example of that...the more samples we add of that rectangular shape in frequency, the closer we get in time to a true Sinc. $\endgroup$ Commented Apr 27 at 11:47
  • $\begingroup$ @MattL. I understand that FIR filter design involves working in both discrete and continuous domains. However, I would like to start with the continuous domain and use the ideal filter responses to gain a better understanding of how to transition to the discrete domain. This will help me connect the theoretical concepts I have studied with the practical aspects of filter creation. But I don't know if the "intuition" I wrote is legit. $\endgroup$
    – minghierid
    Commented Apr 27 at 15:26
  • 1
    $\begingroup$ @minghierid: If there are still some questions left after reading all the answers, you could formulate them by editing your original question. $\endgroup$
    – Matt L.
    Commented Apr 28 at 11:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.