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$
  • 1
    $\begingroup$ it’s never a good idea to set bins to zero $\endgroup$ – Stanley Pawlukiewicz Oct 18 at 1:23
  • 1
    $\begingroup$ dsp.stackexchange.com/questions/6220/… $\endgroup$ – Stanley Pawlukiewicz Oct 18 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$ – M. Farooq Oct 18 at 1:57
  • 2
    $\begingroup$ @M.Farooq proper linear filtering. With a filter! $\endgroup$ – Marcus Müller Oct 18 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$ – M. Farooq Oct 18 at 13:12
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$
  • 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$ – M. Farooq Oct 18 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$ – Laurent Duval Oct 18 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$ – M. Farooq Oct 18 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$ – Laurent Duval Oct 18 at 16:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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