0
$\begingroup$

I am trying to implement fast convolution between a signal and complex Morlet wavelets. To make the result equivalent to a linear convolution, I let the MATLAB fft function zero-pad both the signal and the wavelets up to a length of S + W - 1, where S is the length of the signal and W is the length of the wavelet (kernel).

Here's code used to compute the transform:

% Initialize parameters of signal
signal_length = length(signal);

% Ensure linear convolution
nConv = signal_length + wavelet_length - 1; % Result of convolution S + W - 1
halfK = floor(wavelet_length/2); % Half of the wavelet

%FFT of signal
sigX = fft(signal,nConv);

% Convolution per frequency
for fi = 1:nFrex
    
    % FFTs of wavelets
    waveX = fft(wavelets(fi,:),nConv)';
    
    % Normalize Wavelet
    waveX = waveX./max(waveX);
    
    % Element-wise multiplication
    convres(fi,:) = ifft( sigX.*waveX);
end

% Get powers
powers = abs(convres).^2;

I have attached a picture of the signal (1 second chirp embedded in 5.5 seconds of noise) and a plot of the spectrogram using the power extracted from the first S points of the result of convolution between the signal and wavelets.

Raw signal

Spectrogram

As seen by the images, the power is expected to begin increasing at time = 0 and continue until time = 1 second, but it begins to rise much sooner than expected.

With some consultation, I have managed to recover the correct time location of the signal, but it's not clear to me why it worked this way or if it is reproducible. The thing I changed is instead of letting fft do the zero-padding for the wavelet, I explicitly added zeros evenly to both sides of the wavelet then used ifftshift to move the center of the wavelet to the 1st sample and have the other half wrap around to the end. If I try to remove floor(W/2) points from both ends of the output, the correct time location is un-recoverable. The code example for this process is below:

% Zero pad wavelet
nConv = signal_length + wavelet_length - 1; % Result of convolution S + W -1 
sig_length_diff = nConv - wavelet_length; % Number of zeros to add to wavelet
padlength = sig_length_diff/2; % # of zeros to add on either side of wavelet
zeropad = zeros(1,padlength); % Make zero vector
zeropad = repmat(zeropad,nFrex,1); % Match dimensions for the wavelet matrix (all nFrex = 80 wavelets)
wavelets = [zeropad wavelets zeropad]; % Add the zeros
wavelets = ifftshift(wavelets,2); % Move center of wavelet to the first sample

Then instead of removing half the length of the wavelet worth of samples from either side of the result, I removed all W points from the end. This is not how I understand convolution to add points to the signal though. I would expect half the kernel worth of points to be added to both sides. Why was the entire length of the kernel added to the end?

$\endgroup$
1
  • $\begingroup$ When you normalize the wavelet spectrum, max(waveX) returns a complex value, you should use max(abs(waveX)) instead. $\endgroup$
    – ZR Han
    Commented Nov 15, 2021 at 5:48

1 Answer 1

0
$\begingroup$

There is nothing wrong with frequency domain convolution per se. You are being a bit sloppy with complex numbers. See ZR Han comment

Also look at at:

waveX = fft(wavelets(fi,:),nConv)';

It looks like you are time flipping your zero-padded wavelet. Matlab's transpose operator also performs a conjugate complex operation. It applied to a spectrum that's equivalent to time flipping the time domain signal. Try instead

waveX = fft(wavelets(fi,:),nConv).'; 

Note the dot in front of the transpose operator (which means transpose-only, no conjugates)

See https://www.mathworks.com/help/matlab/ref/ctranspose.html versus https://www.mathworks.com/help/matlab/ref/transpose.html

$\endgroup$

Your Answer

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

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