How to fix a shift issue after taking hamonic mean?

Hi everyone,

In the figure below, the red curve (Mdry) is the harmonic average of blue curve Md by using moving mean of 4 values. Due to harmic mean the red peak has been shifted a bit as shown in figure. How can I fix it? I have attached the data of both Md and Mdry. I am using following code to create a haromic mean (in such a way that the length of both Md and Mdry should remain same -- at the end).

Mdry = zeros(size(Md));
for k=1:length(Md)-3
    Mdry(k)=harmmean(Md(k:k+3));
end
k = length(Md)-2;
Mdry(k) = harmmean(Md(k:k+2));
k = length(Md)-1;
Mdry(k) = harmmean(Md(k:k+1));
k = length(Md);
Mdry(k) = Md(k);

Answer 1

[Mathieu NOE]

hello

is there a reason why the harmonic averaging is needed

here I compared it with a standard smoothing function with gaussain window (odd length) - yellow trace

seems to get similar results in terms of smoothing while better keeping peaks positions

Mdry = zeros(size(Md));
for k=1:length(Md)-3
    Mdry(k)=my_harmmean(Md(k:k+3));
end
k = length(Md)-2;
Mdry(k) = my_harmmean(Md(k:k+2));
k = length(Md)-1;
Mdry(k) = my_harmmean(Md(k:k+1));
k = length(Md);
Mdry(k) = Md(k);

% my suggestion
Mdry2 = smoothdata(Md,'gaussian',7); % NB : odd window length to keep peak position closest to original position

plot(Md)
hold on 
plot(Mdry)
plot(Mdry2)

hold off


function out = my_harmmean(in)

out = numel(in)./(sum(1./in));

end

[Mathieu NOE]

I was a bit slow to simply see the way you did the averaging was indeed creating a shift of 1.5 sample as the k th output is the average of input samples from k to k+3

to avoid this shift , you need to change your code so you pick k samples before and k samples after the running index . So your averaging window has 2*k +1 samples (thus is always odd) and the shift is gone !

here k = 3 , thus window length is 7 ;

try this :

Mdry = Md; % don't bother with the first and last k samples which are not averaged here !

% my suggestion #1
% modified harmmean with odd number of samples and symmetric windowing
for k=4:length(Md)-3
    Mdry(k)=my_harmmean(Md(k-3:k+3));
end


% my suggestion #2
Mdry2 = smoothdata(Md,'gaussian',7); % NB : odd window length to keep peak position closest to original position

plot(Md)
hold on 
plot(Mdry)
plot(Mdry2)

hold off


function out = my_harmmean(in)

out = numel(in)./(sum(1./in));

end

[Mathieu NOE]

hello again

problem solved ?

[Nisar Ahmed]

@Mathieu NOE thanks, I have to check it and then respond you... I have been badly involved some other assigmenets. that's why delayed.

Since I got answer after many dats of posting so become busy in other stuff

[Mathieu NOE]

ok - no problem; take your time !

all the best