# Why result of conv and point wise multiplication has different result?

6 views (last 30 days)
nirwana on 23 Aug 2023
Edited: Bruno Luong on 24 Aug 2023
I am comparing conf function and sum over point wise multiplication of two signal and then plot it. I end up with figure that have different axis length. Literature says that total point convolution is Signal+Kernel-1, but when i check dimension of con_auto is 1x4001. Does it should be
length of cmw (4001)+length of chirpTS (6001)-1?
Does anyone can explain to me?
here is modification code that I use
clc,close; clear
srate=1000; f=[2 8];
t=0:1/srate:6; n=length(t);
% create chirp signal
chirpTS =sin(2*pi.*linspace(f(1),f(2)*mean(f)/f(2),n).*t);
% create complex Morlet wavelet
% wtime is wavelet time, cmw is complex wavelet waveform
wtime=-2:1/srate:2;
w=2*(4/(2*pi*5))^2;
cmw=exp(1i*2*pi*5.*wtime).* exp((-wtime.^2)/w);
% half of the length of the wavelet
halfwavL = floor(length(wtime)/2);
% zero-pad chir --> zero padding dilakukan agar hasil konvolusi mempunyai panjang yang sama
% dengan signal
% CONVOLUTION
for i=halfwavL+1:length(chirpTS)+halfwavL-1
%at each time point,compute dot product
end
% cut off edges
convres =convres(halfwavL:end-halfwavL-1);
subplot(511),plot(t,chirpTS),title ('chirp signal')
subplot(512), plot(wtime,abs(cmw)), title ('morlet wavelet (abs)')
subplot(513),plot(t,abs(convres)), title ('convolution result')
%CALCULATE CONVOLUTION USING CONV Matlab Function
% 'same' - returns the central part of the convolution that is the same size as A.
con_auto=conv(cmw,chirpTS,'same');
subplot(514), plot(abs(con_auto)),
title ('convolution result using conv function')
% CONVOLUTION VIA FREQ DOMAIN
Lconv=length(t)+length(wtime)-1;
convres2=ifft( fft(chirpTS,Lconv).*fft(cmw,Lconv) );
convres2=convres2(halfwavL:end-halfwavL-1);
%line 29 also can be written as
subplot(515), plot(t,abs(convres2))
title ('convolution result via frq domain')

Bruno Luong on 23 Aug 2023
Edited: Bruno Luong on 24 Aug 2023
You have at least 2 errors in your codes
• The convolution by for-loop needs to flip the kernel (Morlet)
• The MATLAB convolution 'same' dimension applied for the first array, so you need to swap arguments
I fix both and here is the result
clc,close; clear
srate=1000; f=[2 8];
t=0:1/srate:6; n=length(t);
% create chirp signal
chirpTS =sin(2*pi.*linspace(f(1),f(2)*mean(f)/f(2),n).*t);
% create complex Morlet wavelet
% wtime is wavelet time, cmw is complex wavelet waveform
wtime=-2:1/srate:2;
w=2*(4/(2*pi*5))^2;
cmw=exp(1i*2*pi*5.*wtime).* exp((-wtime.^2)/w);
% half of the length of the wavelet
halfwavL = floor(length(wtime)/2);
% zero-pad chir --> zero padding dilakukan agar hasil konvolusi mempunyai panjang yang sama
% dengan signal
% CONVOLUTION
for i=halfwavL+1:length(chirpTS)+halfwavL-1
%at each time point,compute dot product
convres(i)=sum(chirpTSpad(i-halfwavL:i+halfwavL).*flip(cmw)); % BUG here, flip cmw
end
% cut off edges
convres =convres(halfwavL:end-halfwavL-1);
subplot(511), plot(t,chirpTS),title ('chirp signal')
subplot(512), plot(wtime,abs(cmw)), title ('morlet wavelet (abs)')
subplot(513), plot(t,abs(convres)), title ('convolution result')
%CALCULATE CONVOLUTION USING CONV Matlab Function
% 'same' - returns the central part of the convolution that is the same size as A.
con_auto=conv(chirpTS,cmw,'same'); % % BUG here, swap signals
subplot(514), plot(t,abs(con_auto)),
title ('convolution result using conv function')
% CONVOLUTION VIA FREQ DOMAIN
Lconv=length(t)+length(wtime)-1;
convres2=ifft( fft(chirpTS,Lconv).*fft(cmw,Lconv) );
convres2=convres2(halfwavL:end-halfwavL-1);
%line 29 also can be written as
subplot(515), plot(t,abs(convres2))
title ('convolution result via frq domain')

Matt J on 23 Aug 2023
In one case, you are using conv(___'same'). The 'same' flag will throw away some of the resulting points.
In the second case, you are doing fft-based convolution, which implements circular convolution, rather than linear convolution. The relation Signal+Kernel-1 only applies to linear convolution.
Bruno Luong on 23 Aug 2023
Edited: Bruno Luong on 23 Aug 2023
I'm lost. You cite a document that is not MATLAB related, or do I miss something?
Just for a fun of OT discussion?
The question is about the difference of the 4th plot; which is obviously not due to small arithmetic errors among different implementations.
Walter Roberson on 23 Aug 2023
There is nothing anywhere in the documentation that restricts MATLAB to produce bitwise identical results to what naive MATLAB-level code would produce for the operation -- and there is a long history of cases where actual implementations produce different results than naive implementations. Actual implementations are permitted to use whatever hardware instructions are appropriate to the situation -- and considering the nature of the Execution Engine in internally producing machine code for execution, the Execution Engine is permitted to produce different hardware instruction sequences for different CPU models.
At the moment, the only restrictions that I can recall as being documented (sometimes documented only in bug reports) are:
• MATLAB restricts itself to not chaining Extended Precision float operations (such as Intel 80 bit registers) -- that it will store back the results of each operation and load again if appropriate or otherwise inject instructions so that each calculation is as-if on 64 bit registers. This restriction might possibly not apply to external libraries it calls such as LAPACK. I do not recall at the moment whether the implication is that MATLAB refrains from using FMA (Fused Multiply and Add) instructions; I think Simulink has a configuration option about whether FMA can be used by Simulink
• When additions are partitioned across multiple cores, MATLAB restricts itself to adding the partial results in a stable order, rather than as the results are available
But, for example, the Windows recommended hardware is "Any Intel or AMD x86–64 processor with four logical cores and AVX2 instruction set support" and at least as of R2023a, notes that at some point in the future AVX2 will be required -- implying that AVX2 is already being used by MATLAB when available, and thus that generally MATLAB gives itself permission to use efficient instructions if the CPU supports them.