Does anyone know about Fourier transforms (in Matlab)?

[mbonsma]

A high school friend of mine is trying to implement FFT in Matlab - I don't know much about the FFT algorithm or signal processing so I'm hoping someone here has expertise in this area.

Our specific question has to do with removing the DC component from a signal - my understanding is that if you subtract the mean, it should remove the peak at 0 Hz, but we're still seeing a peak. Any tips on what this could be and how to address it?

[linamnt]

Someone from my lab does this...I think! writing here to remind myself to ask her!

On Fri, Mar 10, 2017 at 2:58 PM, Madeleine Bonsma notifications@github.com wrote:

A high school friend of mine is trying to implement FFT in Matlab - I don't know much about the FFT algorithm or signal processing so I'm hoping someone here has expertise in this area.

Our specific question has to do with removing the DC component from a signal - my understanding is that if you subtract the mean, it should remove the peak at 0 Hz, but we're still seeing a peak. Any tips on what this could be and how to address it?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/UofTCoders/studyGroup/issues/257, or mute the thread https://github.com/notifications/unsubscribe-auth/AKQ2vlF0VG9m3MY3GoQZARstU022MEAgks5rkat5gaJpZM4MZzUK .

-- Lina Tran PhD Student | Frankland Lab http://www.franklandlab.com/ Department of Physiology University of Toronto

[2ufotofu2]

If the goal is to have a signal in the time domain with no DC component, the simplest way to do this via FFT is taking a FFT, subtracting the low frequencies*, then performing an inverse FFT. An example of this (albeit in python) is in the second answer here. EDIT: I realized this example is a bit misleading, as it subtracts frequencies that have an amplitude less than a threshold, rather than subtracting frequencies below a cutoff frequency.

*Generally, subtracting only 1 frequency (0 Hz in this case) using the above method could cause artifacts when you transform back to the time domain, as discussed here. This might not be an issue for the specific case of 0 Hz, however.

I'm not sure why subtracting the mean doesn't remove the peak at 0 Hz. To try to understand this, try comparing a plot of 1) this mean-subtracted signal to 2) the same mean-subtracted signal, after taking a FFT, removing the 0 Hz, then doing an inverse FFT back to the time domain. If these signals look the same, then I would guess the 0 Hz peak for the mean-subtracted signal is some sort of artifact.

[2ufotofu2]

Here are some (very ugly) plots using that code, using a sine wave shifted up by 1. Note the peaks at 0 Hz and 1 Hz. The sawtooth-looking features in the FFT are some sort of sampling artifact that I can't remember the source of at the moment (probably not relevant to the current discussion):

Compare that with a regular sine wave, which is effectively sin(2\pi x)+1 - mean( sin (2\pi x)+1)--i.e. subtracting the mean removes the 0 Hz component.

This is, of course, basically the simplest example possible, so something more complicated may be going on in your friend's situation. I can't think of what this might be, however. Could you share the data so that we could try to reproduce the problem?

[jladan]

tl;dr: it's probably caused by a low-frequency (slow-changing) part of the signal, like a slow decay.

My masters degree was actually in time-frequency analysis, but it's been a few years, and I'd maybe need to play around with it a little bit to get to the heart of the problem. For the sake of notation, I'll use f for the signal and F = fft(f) for the transformed signal (also, zero-indexing).

First, subtracting the average value definitely does make F(0) = 0. It has no effect on any other component. If you subtracted the average and (in matlab) F(1) ~= 0, then there's a problem in the implementation of either the averaging or the FFT, and there's some debugging to do.

My guess is that there's a very low frequency component in the signal (below the lowest resolvable frequency). Because the DC component is only F(0), a wider peak will remain except for that 0th component. A small linear term would probably do this, like if f(x) = sin(x) + .01*x, with the actual effect dependent on the resolution. This is partly why DC coupling in oscilloscopes are actually implemented with a low-pass filter.

Next, to clarify cgranston's comments.

• the algorithm in the stackoverflow link is a "hard threshold" denoising algorithm. The "soft threshold" algorithm is better for denoising (actually, optimal using the min/max metric and gaussian white noise).
• Removing a single component from F is exactly the same as subtracting a sine wave, so not really an artifact per se, especially for F(0) which just subtracts the average. It's using a rectangular window that causes "ringing" artifacts (exactly what the single-slit experiment demonstrates).
• The jaggedness and wide peak of the sin(x)+1 is due to aliasing. The mesh is set up to be (2pi+pi/50)-periodic. A better mesh when using sin(x) would be x = np.linspace(0, 2*pi - N/(2*pi), N) or x = np.arange(N) * 2*pi/N.

[2ufotofu2]

Thanks for clarifying my comments! All my FFT knowledge is self-taught, so it's nice getting clarification from someone who actually knows what they're talking about!

the algorithm in the stackoverflow link is a "hard threshold" denoising algorithm. The "soft threshold" algorithm is better for denoising (actually, optimal using the min/max metric and gaussian white noise).

How would one implement a soft threshold algorithm? Is it just as simple as having a gradual tapering off of frequency amplitudes below f_0, instead of having a hard cutoff at f_0?

Removing a single component from F is exactly the same as subtracting a sine wave, so not really an artifact per se, especially for F(0) which just subtracts the average. It's using a rectangular window that causes "ringing" artifacts (exactly what the single-slit experiment demonstrates).

By "artifact" I was referring indeed referring to ringing (throughout my post I used "artifact" as a perhaps not correct catch-all term for anything undesirable).

If a single component is removed from F and nothing else is done, isn't this effectively a rectangular window? From the dsp stackexchange post I quoted:

...Thus zeroing a bin will produce the same result as subtracting that sine wave, or, equivalently, adding a sine wave of an exact FFT bin center frequency but with the opposite phase. Note that if the frequency of some content in the time domain is not purely integer periodic in the FFT width, then trying to cancel it by adding the inverse of an exactly integer periodic sine wave will produce, not silence, but something that looks more like a "beat" note (AM modulated sine wave of a different frequency).

The jaggedness and wide peak of the sin(x)+1 is due to aliasing. The mesh is set up to be (2pi+pi/50)-periodic. A better mesh when using sin(x) would be x = np.linspace(0, 2pi - N/(2pi), N) or x = np.arange(N) 2pi/N.

I'm a little confused on this. To avoid aliasing, is the idea to have the mesh have the same periodicity as the signal? Also, how exactly is this aliasing? My understanding of aliasing is that frequencies in your signal that are less than half of your sampling frequency appear as lower frequencies.

Changing the mesh to your second suggestion does make the FFT look much better, thanks!

[jladan]

How would one implement a soft threshold algorithm? Is it just as simple as having a gradual tapering off of frequency amplitudes below f_0, instead of having a hard cutoff at f_0?

Basically, with soft-thresholding, you perform a hard threshold, but also subtract the threshold value from the magnitude of all components. Donoho wrote a paper about it in 94 using orthogonal wavelets (which are probably better for the purpose).

Also, how exactly is this aliasing? My understanding of aliasing is that frequencies in your signal that are [more] than half of your sampling frequency appear as lower frequencies.

Yes, frequencies that are higher than the Nyquist frequency are not distinguishable from the frequencies that make up the basis (0, ... , N/2). That means they are picked up as an additional value on top of the basis frequencies. That phenomenon is aliasing -- the non-basis frequencies respond to the basis frequencies' name.

As a result, frequencies below the Nyquist frequency also exhibit aliasing if they are not perfectly in the basis. One way to think of it is that they respond a little bit to all the closest frequencies. Another is to notice that the full Fourier series has infinite non-zero coefficients, because $ a_k = \int_0^{2\pi} \sin(fx) \sin(kx) dx $ is not zero if $f$ is not an integer.

[jladan]

@mbonsma has your friend figured out their problem yet?

[2ufotofu2]

@mbonsma Feel free to interrupt my thread-derailing!

Yes, frequencies that are higher than the Nyquist frequency are not distinguishable from the frequencies that make up the basis (0, ... , N/2).

What is between 0 and N/2, i.e. what determines the basis frequencies?

I haven't heard of basis frequencies before, only basis functions (e.g. sine and cosine are the basis functions for Fourier series). Are basis frequencies related, and can you recommend a source for reading more about them?

[jladan]

I haven't heard of basis frequencies before, only basis functions.

Sorry, I was just using my own kind of shorthand, and skipped writing 2pi. When I said basis frequencies, I was referring to the frequencies of the basis functions. The discrete Fourier transform using the basis functions exp(i2pi/N kx) for k in (-N/2-1, .., N/2) or (0, .., N-1). I was referring to the value of k, or the index of the basis function.

[2ufotofu2]

Sorry, I was just using my own kind of shorthand, and skipped writing 2pi. When I said basis frequencies, I was referring to the frequencies of the basis functions. The discrete Fourier transform using the basis functions exp(i2pi/N kx) for k in (-N/2-1, .., N/2) or (0, .., N-1). I was referring to the value of k, or the index of the basis function.

Got it, thanks for explaining.