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.
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?
max(waveX)
returns a complex value, you should usemax(abs(waveX))
instead. $\endgroup$