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

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
chirpTSpad = [zeros(1,halfwavL) chirpTS	zeros(1,halfwavL)];

% CONVOLUTION
convres	=zeros(size(chirpTSpad));
for	i=halfwavL+1:length(chirpTS)+halfwavL-1
		%at	each time point,compute	dot	product
		convres(i)=sum(chirpTSpad(i-halfwavL:i+halfwavL).*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(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')

Answer 1

[Bruno Luong]

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
chirpTSpad = [zeros(1,halfwavL) chirpTS	zeros(1,halfwavL)];

% CONVOLUTION
convres	=zeros(size(chirpTSpad));
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')

Answer 2

[Matt J]

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.

[Walter Roberson]

conv() is permitted to decompose the convolution, and to use fft() based methods, and to use hardware SIMD (Single Instruction Multiple Data) instructions, and to partition the problem over multiple cores. It was never guaranteed to return exactly the same results as naive looping.

[Bruno Luong]

@Matt J "In the second case, you are doing fft-based convolution, which implements circular convolution,"

Incorrect, since @nirwana uses fft with 0-padding to input arrays so it is circular but on the bigger arrays and the circular effect is simply NOT there.

[Bruno Luong]

@Walter Roberson "... use hardware SIMD (Single Instruction Multiple Data) instructions ..."

Where is that in the code?

[Walter Roberson]

@Bruno Luong

conv() mostly involves conv2() to do the work.

conv2() is built-in so we cannot see the code for it. The internal code could be highly optimized. For example it might call Intel® oneAPI Math Kernel Library VS routines for convolution https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-0/convolution-and-correlation.html

[Bruno Luong]

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]

@Bruno Luong

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.