Convolution via frequency domain multiplication
Вставка
- Опубліковано 6 вер 2024
- Is time-domain convolution too slow? (Yes it is.) Learn how to do lightning-fast convolution in the frequency domain. This will also help you understand that wavelet convolution is really just filtering.
The video uses files you can download from github.com/mik...
For more online courses about programming, data analysis, linear algebra, and statistics, see sincxpress.com/
Honestly a great tutorial. I do radar signal processing and this is a wonderful refresher for pulse compression techniques.
Awesome, I'm glad it was helpful :)
Wow, this is really enlighting. Have to watch it on the PC again an I will try to make a Java implementation.
I spent way longer than I want to admit implementing this in javascript, and learning the lesson that multiply in this context refers to multiplication of complex numbers and not just regular old multiplication lol
Yes, convolution is tricky when you first learn about it...
Amazing lecture!
if you apply the wavelet transform to the audio signal of the video, you can actually hear the spinning cycle of the washing machine :)
Ha, yeah that's possible. I have another recording where you can hear a helicopter outside! I do these recordings in my apartment (which is in the middle of the city), so there is sometimes "real world" background activity.
@@mikexcohen1 a meta-application of your content, I love it. Great lectures by the way; have helped me understand the wavelet transform in a clear and concise way.
wow it's amazing how much faster this approach is. I was messing around with convolutional reverb in matlab and using just a 10sec audio file and the conv command would have my pc locked down for a min as it was calculating the convolution between the impulse response and input sound file. But doing it this way took literally 3 secs if that lmao.
lol just did the reverb on a 3min audio file which using the conv command would just crash my matlab but it was able to do it so fast.
Nice ;)
Great lecture. Do you have any videos on ERS/ERD? Also on Classification using machine learning? I am always confused regarding the choice of features? Use of common spatial patterns etc. It will be helpful for neural signal processing learners if you include those in your book and upload video lectures.
Hi Debanjan. ERD/S are just names for relative power de/increases. The videos on extracting power from Fourier and wavelet coefficients, and the videos on power normalization and choice of baseline, should be helpful. Unfortunately, I don't have any videos directly on classification or machine learning. However, many linear multivariate classification algorithms rely on generalized eigendecomposition, and I have a few videos about that.
Mike
nice explanation
Hi Mike, thank you for the superb video. I have a question:
When you plot the power spectrum of the data, you have:
2*abs(dataX(1:length(hz))/length(data));
You mentioned the factor of 2 is included because of the positive and negative frequencies. I was wondering what is the divisor of "length(data)" for? Also, when you plot the power spectrum of the wavelet, there is only:
abs(cmwX(1:length(hz));
Why are the factor of 2 and the divisor not included?
Another thing: when you use line 89 to cue 1/2 of the length of the wavelet from the beginning and from the end, it looks like the IFFT result starts and ends one data point earlier than the intended result. I think the code for line 89 should be:
conv_res_timedomain = conv_res_timedomain(half_wav:end-(half_wav-1));
Hi Chao. I have a video about extracting correct units from the result of FFT. Basically, there are two normalizations you need: divide by N and multiply by 2. The former is due to the summing in the Fourier transform, and the latter is due to including the positive and negative frequencies (which are the same for a real-valued signal, hence multiply by 2).
Thank you! Now I finally understand that.
Great job!
Quick question, suppose I have a signal size 16 and a kernel size 3. How should I pad the kernel? At the end? Or some padding in the begin and the rest at the end, so that the kernel is centralized?
You don't need to zero-pad the kernel yourself; you would let the fft function do this for you. Thus, something like fft(kernel,N_conv), where N_conv is the length of convolution (N+M-1).
Thanks again for the great lecturelets. Is it possible to extract the standard deviation of a frequency domain Gaussian from it's time domain Morlet parameters?
bw,
Eivind
Best to estimate it empirically. See this paper: www.biorxiv.org/content/biorxiv/early/2018/08/21/397182.full.pdf
Thanks! According to your paper I can derive the Morlet wavelet from a frequency domain Gaussian. Simple, yet brilliant. As a side note, it seems to be a strong "sd=log2(frequency)" relation across frequency to time domain.
Indeed it is a logarithmic relationship. However, I still recommend empirically measuring it, because sampling rates and data lengths will impact the actual filter that is applied to the data, and I think it's more important to know the actual characteristics than the analytically specified parameters. You can tell that I'm more an empiricist than a theorist ;)
Great, there are some minor "inaccuracies" in my results that I can't explain, and they are probably related to your point. I've been trying to measure band-power for varying frequencies 'x' using 1) standard FFT , 2) constructing a Morlet with central frequency 'x', 3) constructing time domain Gaussian with specified SD. When I set 2^(log2(freq)-3) as SD in frequency domain Gaussian or as bandwidth in standard FFT, I get almost exactly the same result as using a Morlet, across all measured frequencies. Exciting application, but I need to fail-proof the math.
Hi, I have a question about the amplitude normalisation of the wavelet in the frequency domain. When I normalise using:
cmwX = cmwX ./ abs(max(cmwX))
the amplitude of the convolution result is half what I would expect (eg input = sine wave of amplitude 1, output = real part of convolution result with amplitude of 0.5)
Do you have any idea why the result is halved?
(NB I used abs(max(cmwX)) instead of max(cmwX) because of weird phase shifts in the result when normalising with the latter, but the value of the max is x + 0i so I don't see how this would affect the amplitude...)
ps thank you so much for these videos! Working through these and your book has been immensely helpful.
Hi Soraya. Great that you are testing the methods yourself! That's the best way to learn. Briefly: this happens because the amplitudes in the Fourier transform get split between the positive and the negative frequencies. So the "missing" half of the amplitude is in the negative frequency spectrum. A full recovery of the amplitude spectrum requires combining the positive and negative frequency coefficients, which in practice means you can double the positive-frequency amplitudes. You can learn more about that in this video: ua-cam.com/video/Ls7AvuZG4kI/v-deo.html&ab_channel=MikeXCohen
And I have much deeper explanations about that in this course: www.udemy.com/fourier-transform-mxc/?couponCode=MXC-FOURIER10
Hope that helps!
Mike
Hi Mike, that's exactly what I was looking for - thank you! And thanks for the speedy response - this factor of 2 has been vexing me for a little while now.
I have one more question about wavelets: in Analyzing Neural Time Series Data (chap 13) you define a frequency specific scaling factor:
A = 1/ (s*sqrt(pi))^1/2 [where s = gaussian std dev]
when creating the complex Morlet wavelets.
The cmwX = cmwX ./ abs(max(cmwX)) normalisation has been necessary to get correct responses out of test sine waves, but removes A from the convolution (in my hands anyway...). Is this A scaling specifically for neural data to tackle the 1/f noise?
Thanks again for you help!
Soraya
The short answer is that normalizing the wavelet in the frequency domain is much better and is always accurate, while normalizing in the time domain is a huge pain and most people do it wrong. You can read a longer answer in this thread: sccn.ucsd.edu/pipermail/eeglablist/2016/011536.html
I didn't describe this well enough in the book, which I regret. A revised edition (if there is one!) will definitely cover this in more depth.
Ah I see- thanks again!
I like your video but you left out something very important. You forgot to tell people that you have to TIME REVERSE one of the functions in the time domain before you do the convolution. I do have a question though. When we do the convolution in the frequency domain, do we have to time reverse one of the functions in the time domain before we take its fourier transform?
Hi b r. In fact, the kernel is reversed only in time-domain convolution. When you implement convolution via spectral multiplication, you don't reverse the kernel. And if you do reverse the kernel before taking its FFT, then that result is equivalent to cross-correlation in the time domain.
Or you can only ever use symmetric kernels; then you don't need to worry about time-flipping ;)
Hi Mike, quick question, lets say you were interested in graphing the phase angle with respect to time, as you demonstrated in this video. Are thee any added benefits of using a complex morels wavelet to obtain this information as opposed to the normal complex sin wave via the FFT? If your signal is non-stationary, will the morels wavelet give you any type of added benefit because its better adapt at dealing with such signals? Thanks.
Hi Brit. If you use a single Fourier coefficient, then you'll just see a pure phase angle time series. If the signal has some fluctuations (non-stationarities), then you'll need more than just a single Fourier coefficient, because the nonstationarities are represented in the side-lobes, not in the spectral peak. I talk about this in the video on stationarity and effects of violations. A Morlet wavelet is Gaussian in the frequency domain, and thus captures some of the side-lobes, meaning it would also show some of the temporal fluctuations in the phase. Morlet wavelets are really great for brain electrical signals, but I don't have enough experience in other areas to say that they would necessarily be great for every type of signal. Nonetheless, it's worth trying.
Mike
Thanks for getting back to me as always Mike, I'll take a few more looks at that particular video and keep tinkering away at it. Would you mind if I emailed you with my findings relating to this particular question once I get some results?
You're always free to email me, particularly if it's about data analysis ;)
Mike, thank you so much for your superbly clear lectures. One question, how does one pad out the odd-length wavelet to an even-length power-of-two ready for the FFT? Surely, convolution via frequency domain multiplication must force the wavelet to be even in length. Thanks so much.
Hi Julian. Why does convolution need to have an even length? As long as the length is appropriate (length of kernel + length of signal - 1), it will work. Having a length of 2^n is not necessary for modern implementations of the FFT.
Hi Mike, thanks for your reply. I gathered in order to do a time-domain convolution via the frequency domain, the two frequency representations (and therefore time) need to be the same length so that we can meaningfully dot-multiply their frequency spectra before transforming back to the time domain. I've therefore zero-padded my otherwise shorter and odd-lengthed kernel out to the same length as the signal (which also respects the original signal+kernel length -1). I'm using an FFT that requires a 2^N size; I'm not using Matlab. I'm doing that for speed, but I guess I could consider alternative DFT transformers that don't have this requirement if necessary. I do though see that in your video, your nConv is 1664 long, which is indeed an even number, so I guess Matlab has to make a decision as to how to automatically pad your odd-lengthed kernel array out to an even length. I'm guessing, aside from doing an odd-length DFT, Matlab is padding so that the original central sample is now as close to the array centre as possible. Talking of which, I also had to shift the entire kernel by one-half its length (like cutting it in half and swapping the left and right halves) in order to get a Gaussian all-positive frequency domain signal. If I don't do this first, my Gaussian is correct except every other frequency value (Re) is negative. Thanks again for such a wonderful series of lectures and also for taking the time to help with my query/queries; it's very much appreciated. Julian
Ah, I see what you mean. In that case, you should compute the length of convolution (M+N-1), then find the next-highest 2^k for the FFT. Then after the IFFT, you trim the extra zeros off the end, and then you clip off the "wings" of convolution (1/2 the length of the wavelet from the beginning and end). Don't actually implement a loop-based DFT unless it's for really small time series. I talk about this 2^k stuff in my ANTS book, and I think it's also in the MATLAB code that accompanies the book, which you can download from sincxpress.com.
Great, that's exactly what I'm doing, but my query was more about the positioning of the odd-lengthed-and-centered kernel within the longer even-lengthed array. I see now, after some further research, that rather than it needing to be centred, it should be shifted and wrapped left such that the array starts with the center (t=0) value, and then filled with the rest of the second half of the kernel, and then the first half in order for the frequency domain transform to be all positive values of the Gaussian. I think Matlab is doing this shift automatically like it pads with zeros. Again, thanks so much. Your lectures are really superb and thank you for the link and extra information. Julian
Instead of triming the 2 vectors you want to compute, you can 0-pad them both to make them the same length. The 0 pad you perform should be before applying fft to the 2 signals. Then you can multiply them and perform ifft and that's it
Can this technique be used for spectral analysis? (spectrograms / CQT spectrograms). From what I can see, this is a fast way to perform bandpass filtering - not time-frequency analysis
That depends on how you set it up. The real part of the result of convolution with one wavelet gives you a narrowband-filtered signal. If you use a complex-valued wavelet you can extract time-varying power. And then if you repeat for many wavelets over a range of frequencies, then it will give you a spectrogram.
@@mikexcohen1 thanks a bunch. I've been trying to figure out how to extract the time-varying power. I was doing inverse real transforms, which could only result in a varying sine wave
@@mikexcohen1 From my experiments, both components of ifft(fft(morlet_wavelet) * fft(signal)) still oscillate (for a constant sine wave signal). I was expecting a constant stream of 1s (representing constant volume of 1 for this frequency)
You need to make sure to set up the N's of the fft correctly, accounting for zero-padding. Ringing effects can come from a variety of sources, including poorly specified wavelets or spikes in the signal. You should check out my time-frequency playlist in the ANTS2 module.
@@mikexcohen1 Thank you very much. That ANTS playlist answered SOOO many questions for me. I was able to achieve the volume (stream of 1s for the sine tone) using the Hilbert transform
Going from time domain to freq. domain results in a complex spectrum but you only show e.g. the magnitude or real part in the figure.
In MATLAB, the fft result is an unscaled, double-sided complex spectrum which you leave as is (as opposed to making it single-sided, properly scaled with DC and Nyquist comps. fixed); this is ok but, again, it's inconsistent with what you are effectively showing in the figs.. Would be nice also if you didn't kind of gloss over the normalisation part of the kernel FFT; maintaining the correct scaling in the output is very important in many (engineering) tasks. Thanks
Hi Ben. That's right -- the multiplication is between the full complex-valued spectra. Showing only the power spectra here is done to facilitate intuition (visualizing a complex-valued Fourier spectrum is difficult). Normalization is indeed important, and I discuss it in different videos. But normalization can be tricky, because different normalizations can be applied, depending on the goal.
@@mikexcohen1 Hi. Thanks for replying. We're trying to do something different which is to perform freq. dom. multiplication by time domain conv.. The challenge is, this, of course, means a kernel of the same length as the orig. time data file. Since we're engineering, we also need to maintain the units. The spectrum that the kernel comes from is am FRF of nominally 1 (but there are diffs. to 1 at certain freqs. which is the reason we are doing this) with zero phase. On the scaling part, I think I've worked out that we should divide thro' by the integral of the FRF (which, if 1, is basically the length of the vector) - sound right? Re. getting a result of the same length as the input file, MATLAB's conv 'full' seems to suggest we should take the front half, not the middle (which you would get from 'same') but I'm not sure why yet. Might need to implement convolution integral in my own code. Would be great if you had any further thoughts. Thanks, Ben
Yes, usually the normalization is either based on the integral of the kernel or the peak of the kernel (forcing either to be 1). Regarding the convolution details, I think you'd still want to pad both sides and then take the middle. It might depend on whether you have Nyquist or DC in the middle of the spectrum. I'm actually not really sure on this point, but my advice is to run it in a small simulation where you know what the ground truth needs to be after convolution.