I'm performing frequency convolution by multiplying FFT coefficients by values computed from a Laplace equation (e.g. H(s)=wc/(s+wc)). The pseudo-code looks like:

X=fft(data);

for ii=1:N/2  % loop over positive frequency bins up to, but not including, Nyquist

    s=1i*2*pi*freq(ii);    

    H(ii)= ...some s-domain equation...

    Y(ii)=H(ii)*X(ii);    

end

I complete the negative-frequency half of the spectrum by mirroring the conjugate of Y(ii) starting at index N/2+2. When I complete Y, then I take an IFFT(Y) to get a time series.

My question is, how to compute the filtered value at the Nyquist bin... that is, Y(N/2+1)? The H(s) value will be a complex number. The original time series data signal is real, as is the expected IFFT time series. 

So far I've been computing it as:

    Y(N/2+1)=X(N/2+1) * abs(H(N/2+1)), 

    where abs(H(N/2+1))=sqrt( (real(H(N/2+1)))^2 + (imag(H(N/2+1)))^2 ).

but, perhaps it should just be

    Y(N/2+1)=X(N/2+1) * H(N/2+1)

then take the real part of the time series after IFFT(Y). Or, something else? I'm not sure. Does anyone know the correct way to filter at the Nyquist bin?

The very first thing to remember here is that you are doing some approximation of some real-world continuous filtering.

The second thing to remember is that if your \(H \left ( s \right ) \) has a significant amount of response left at Nyquist, your approximation of whatever real-world thing you're doing isn't very good.

Put these two things together, and it means that you are either sampling too slowly for your ambitions, or you can do whatever you want at Nyquist.  I vote for just using the real part of \( H \left ( j\omega \right ) \) at Nyquist, or, preferably, verifying that \( H \left ( j\omega \right ) \) is less than 0.01 or so, and zeroing it out.

Thanks Tim.

The sampling rate is fixed by the system, and has a dominant tone centered at Nyquist. The magnitude of this tone is often very much higher than all other tones in the spectrum. I expect it should be centered exactly at the Nyquist bin (or pretty darn close to it). So even if H(s) is much lower than 0.01, H(s)*X(s), the filtered spectrum, is often dominated by the tone at Nyquist, and it's important to represent it in the time domain, as least within the limits of the algorithm. 

Also, the filter parameters (cutoff freqs), are system dependent, and are independent of choice of sampling frequency. This sounds like it's a little different than the systems you work with.

The Nyquist tone, as I call it, includes the Nyquist bin plus a few bins on either side of it that form an apron. As an example, if I list the FFT coefficients, X(i), and filter coefficients, H(i) for the positive-frequency apron plus Nyquist bin, I have: 

H(i):

      6.480334208022409e-04 - 2.544825089250270e-02i

      6.480332663991106e-04 - 2.544824786276553e-02i

      6.480331119960354e-04 - 2.544824483302910e-02i

      6.480329575930155e-04 - 2.544824180329338e-02i

      6.480328031900510e-04 - 2.544823877355839e-02i

X(i):

     -1.328191551980880e-09 - 8.160753005573973e-09i

     -1.319711956620663e-07 + 5.187901186851907e-08i

      1.598673069238291e-06 - 2.605737644498312e-07i

     -5.464408831672277e-06 + 4.206131308438971e-07i

      8.006485160682597e-06 + 0.000000000000000e+00i

listed in order for the following bins:

      bin N/2-3

      bin N/2-2

      bin N/2-1

      bin N/2 

      bin N/2+1 (Nyquist)

As a baseline, if I do NOT apply H(s), and simply IFFT the X coefficients, everything looks fine in the time domain (the signal's amplitude bounces up and down at the same peak-peak values). But after I apply the H(s) filter then IFFT the result, I get a time series with a low frequency envelope. I don't know why.

"The sampling rate is fixed by the system, which uses the rising and falling edge of a clock signal. It cannot be changed. If the clock signal is 100 MHz, the sampling rate is 200 MHz."

How square is the clock signal? If it goes "lub dub .. lub dub like my heart, that will introduce artifacts due to the phase modulation. Spectral purity gets to be important.

So the really high frequency content is important to you, and you're low-pass filtering the data to get rid of the high frequency content?

Why?

I'm feeling like I'm the trickee of a trick question.

Yes. I understand your perspective. Note that we're analyzing noise in the system. In the real world, it can appear at any frequency and any magnitude. A system designed with a low-pass filter, for example, can still have noise beyond its cutoff frequency. Also, the filter's roll-off isn't brick wall, but rather like -20 dB/dec. So even one decade beyond the cutoff, if noise is present with a sufficient magnitude, it can effect system performance. My goal is to model the system, and the noise, to quantify that effect.  

Use the real part then, or try various things.  Nothing is going to be quite right up there, thanks to your choice of how you represent things.

Thanks Tim, Could you help me pinpoint where the algorithm breaks down for Nyquist? (the more specific the better) Any lower-frequency tone appears to work just fine (for example, a tone centered a dozen bins lower than Nyquist). What assumption am I breaking at Nyquist?

I agree with Tim - there should be no amplitude at the Nyquist limit when simulating an analog system.  As you have discovered, it is impossible to specify both the magnitude and phase at that point.  It's the first "out of range" frequency for the sampling theorem, being both +fs/2 and -fs/2.  In other words, fs/2 and -fs/2 have been aliased on top of each other, canceling the imaginary part to zero for any real signal.

What I would do is either (1) multiply the existing system H(s) by an analog lowpass filter G(s) that has a null at fs/2, making sure that the FFT length is longer than h * g (to avoid being undersampled in the frequency domain), or (2) apply the bilinear transform to H(s) and accept the frequency warping that comes with that.  I normally do (2) with an oversampling factor sufficient to make the frequency warping undetectable in my application.

Thanks JOS, 

Since H(s) isn't real at the Nyquist bin, I wasn't sure whether the filtered FFT coefficient at that bin should be real or imaginary, before IFFT'ing the spectrum. Thus, the post. Your approaches to improve the processing I'm sure are better than what I'm doing. However, I'm still trying to understand the limits of the process I'm using, before making it more complex, as I'm not that much an expert in filtering. That said, could you explain a little more the statement, "there should be no amplitude at the Nyquist limit when simulating an analog system"? Also, could you clarify, is that no amplitude in H(s) at Nyquist? 

Yes, I mean simply that we need H(j 2 pi f) = 0 for all |f| greater than or equal to fs/2 (half the sampling rate = Nyquist limit).  Otherwise, we are violating the conditions required for the sampling theorem to hold.  Since information is lost at fs/2, there is no way to invert the sampling operation and recover the original continuous-time system.

Thank you JOS,

When I perform an IFFT of the Nyquist tone (including Nyquist bin and its related apron), it reconstructs the original signal fine. When I apply H(s) and repeat, it doesn't. Are you saying that using H(s) introduces an additional requirement at Fs/2 (beyond the normal requirements for FFT, IFFT) to reconstruct a signal? If so, could you explain in more detail? If there is no signal above Fs/2, then there is no consequence of having H(s) be non-zero above Fs/2 (in practice, H(s) above Fs/2 will never be used in my algorithm). But why must H(s) be zero at Fs/2? Isn't Fs/2 still useful for reconstruction?

We often get sloppy about limits when thinking in terms of continuous variables. Remember, we count frequencies with an FFT, we don't measure them as such. Fs/2 is the lowest frequency that is too high to meet all the sampling criteria. Don't confuse "<" with <=". I forget who said, "Funny things hppen at the margins."

Thanks jya, Can you explain "Fs/2 is the lowest frequency that is too high to meet all the sampling criteria"? 

I thought Nyquist's sampling criteria is that the sampling frequency (fs) must be at least twice the highest frequency (fc) contained in the signal, or, fs >= 2*fc. Isn't the Nyquist bin at fc, and therefore satisfies this criteria (with the equal sign)? 

There is an acquisition time. The closer you get to DC or Fs/2, the longer you have to squint at the signal to understand it. The relation is essentially reciprocal, so that at DC or Fs/2, you need to look forever.

It's easier to illustrate the DC case with words, but whatever is true of f is also true of Fs/2-f. Here's an illustrative problem:

Fs = 20 KHz,  f = .1Hz

How many samples are needed to reliably tell that f is not DC? Ten times as many would be needed for .01 Hz.

So, is it that the Nyquist bin (e.g. Fs/2) satisfies the sampling criteria, but only if the measurement is infinitely long? If so, then does that mean the Nyquist bin contains no useful information for reconstruction (or otherwise, how to interpret the data in its bin)?