1
$\begingroup$

I have a real one dimensional signal s (light absorbance in a flow cell), which has significant noise and some periodic noise after performing a deconvolution of $S$ from $S_o$. Basically fft($S$) was divided by the fft($S_o$). Let us call the output as $R$ vector. In order to remove the periodic noise, I replaced the noise region in $R$ with 0+0i in a given range. The ifft($R$) is complex. When I plot t and ifft($R$), MATLAB warns that the imaginary parts of complex R arguments ignored . The output looks as desired, free from noise.

Forum members mention that this can arise due to numerical precision, however the imaginary part in ifft($R$) is on the order of 0.0006i.

Is there is a better window (simple one) than this rectangular window in the frequency domain?

Even in this simpler example which one can test is a simpler code, generates a complex inverse output.

t = [0:1/80:60]'; % Time
x = sin(2 * pi * 2 * t) + sin(2 * pi * 0.05 * t); % Signal
X = fft(x);
X([1:6, 4780:4801]) = 0; % Removing low frequency 0.05 Hz
z = ifft(X); % **z turns out to be complex as well**
figure(1)
plot(t, z)
hold on
figure(2)
plot (abs(X))

Thanks.

$\endgroup$
5
  • 1
    $\begingroup$ it’s never a good idea to set bins to zero $\endgroup$
    – user28715
    Oct 18, 2019 at 1:23
  • 1
    $\begingroup$ dsp.stackexchange.com/questions/6220/… $\endgroup$
    – user28715
    Oct 18, 2019 at 1:25
  • $\begingroup$ Thanks. I had read that answer. What is a good alternative to remove to undesired frequencies in the frequency domain if zeros is the worst case scenario (rectangular window). $\endgroup$
    – AChem
    Oct 18, 2019 at 1:57
  • 2
    $\begingroup$ @M.Farooq proper linear filtering. With a filter! $\endgroup$ Oct 18, 2019 at 8:12
  • $\begingroup$ @MarcusMüller, I am coming from pure chemistry side, so I never had formal DSP courses. This is all a self-learning process over a year or so arising as per need basis from books and online lectures! $\endgroup$
    – AChem
    Oct 18, 2019 at 13:12

1 Answer 1

3
$\begingroup$

Essentially, your code does not respect the inherent Hermitian symmetry of the output of the FFT. Here, your signal is odd-sized $2K+1$. Hence, this FFT yields a complex vector of coefficients $d$ (real) and $a_k$ (generally complex), arranged as:

$$ \left[d,a_1,a_2,\ldots,a_K,\overline{a_K},\ldots,,\overline{a_2},\overline{a_1} \right]$$

If you want to zero-out some coefficients, you should remove then by conjugate pairs (except for $d$), like in

$$ \left[0,0,a_2,\ldots,a_K,\overline{a_K},\ldots,,\overline{a_2},0 \right].$$

Note: for even sizes, a real coefficient $n_{K+1}$ is inserted between $a_K$ and $\overline{a_K}$.

Your line:

X([1:6, 4780:4801]) = 0

does not satisfy this symmetry. Try instead this one (only valid for odd length):

nACBin = 3;
X([1:nACBin+1, end-nACBin+1:end]) = 0; % 

and you will get a much lower imaginary residue (here due to numerical issues). Here, it somehow works (see below) because of the very simple context (two sines with lucky interval length), but this is rarely recommended in general:

Sine signals and rectangular window Fourier filtering

Especially with short 0 strides, this will convolve your data with a cardinal sine, both broad and wiggling.

As said in comments, you are considering prototyping a filter in the frequency domain using windows. This is a well-covered domain in DSP. However you can avoid those complexities by first looking at filtering in the time domain (Butterworth, Chebyshev, Cauer, Papoulis filters) that already have well-designed implementations. Other filters can be useful for peaks, like polynomial interpolators (Savitzky-Golay, for instance).

Side note: it probably would be more efficient (though a bit more complicated) to embed those frequency constraints into the deconvolution, instead of trying to remove the artifacts afterward.

$\endgroup$
4
  • 1
    $\begingroup$ Thanks for this great insight on reducing the imaginary components. The typical time domain filters like the numerically simulated RC introduce tailing in Gaussian peak shapes. I have published in detail on some common time domain filters used in instrumental chemistry. Moving average also broadens the peaks, that is why I was considering removing those specific frequencies in FT which were causing ripples and spikes. I am coming from pure chemistry side, so I never had formal DSP courses. $\endgroup$
    – AChem
    Oct 18, 2019 at 13:19
  • $\begingroup$ You are welcome. Sure, traditional RC filters broaden peaks. If your data is peak-like, polynomial filters, like Savitzky-Golay, designed by chemists, can be a good start. $\endgroup$ Oct 18, 2019 at 15:27
  • 1
    $\begingroup$ Laurent, it is actually chromatographic data. I believe you will understand the problem. I wanted to remove HPLCs broadening by Fourier deconvolution in the chromatogram. The chromatographic data was deconvoluted by the instrument response which led to a some spurious frequencies. Just two days ago, I was telling students about Savitsky-Golay's history. Golay was a brilliant electrical engineer and perhaps Savitsky was a chemist. SG is not as great as a centered weighted moving average (Hamming etc.) in terms of removing noise in chromatograms. $\endgroup$
    – AChem
    Oct 18, 2019 at 15:32
  • $\begingroup$ There are generalizations of SG. Plus, sparse deconvolution might by adapted if you have ideas about the shape of the broadening $\endgroup$ Oct 18, 2019 at 16:17

Your Answer

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

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