フィルターのクリア

n point symmetric weighted moving average filter

6 ビュー (過去 30 日間)
Luke Pretzie
Luke Pretzie 2016 年 10 月 27 日
コメント済み: Star Strider 2016 年 10 月 31 日
The following is a hard coded 3-point weighted symmetric moving average filter:
ECG(:,1) = time; %initialization of the data, rather crude, I have yet to streamline it: the important part is N.
%Attached is the ECG data: The first column is the time vector, and the second column is the voltage vector of the signal.
ECG(:,2) = Amplitude;
N = length(ECG);
A = ECG(:,2);
time = ECG(:,1);
% 3 point Weighted Symmetric filter
for i = 2 : (N-1) %Filtering starts at the second point since you cannot gather data from point 1 since there is no previous point.
%We also ignore the last point for similar reasons
W3(i-1)= (A(i-1) + 2*A(i) + A(i+1))/4; %Adds weights of one to the first and last points in the filter, and increases the weight value by 1 as you move
%closer to the center point. Thus a 5 point filter would have weights of 1,2,3,2, and 1, respectively.
%You then have to divide the entire thing by the sum of the coefficients, which in this case is 1+2+1 = 4.
%If this was a five point filter, it would be 1+2+3+2+1 = 9 that the expression would be divided by.
end
plot (time (2 : (N-1)), W3, 'g') %Then all that's left is to plot the data, with necessary labels and a title
xlabel('Time (sec)');
ylabel('Amplitude (Volts)');
title('3 Point Weighted Symmetric Filter');
So my assumptions for how an n-point weighted symmetric moving average filter would function are as follows:
% n point Weighted Symmetric filter: Let's assume n is set equal to 3, so that it is easy to follow from the first code
reach = n - ceil(n/2) %How far out the filter "averages". For instance, a 3 point filter would have 3-ceil(3/2) = 3-3 = 1 neighboring point both before
%and after the point of interest that it would be averaging over.
for i = ceil(n/2) : (N-1)
WN(i-reach)= (reach*A(i-reach) + (reach+1)*A(i) + reach*A(i+reach))/((3*reach)+1); %??? Yeah, it gets tricky for me at this point.
end
plot (time (ceil(n-2) : (N-1)), W3, 'g') %Sure enough, our final goal is to plot this data
xlabel('Time (sec)');
ylabel('Amplitude (Volts)');
title('n Point Weighted Symmetric Filter');
%}
My final goal is to create a weighted symmetric moving average filter that has a modular number of points over which it can average. The part that really gets me is the weighting itself, and while I'm sure that a nested for loop of some kind would do the trick, I can't see how I would even start something like that.
Thank you for taking the time to inspect my question, any feedback would be greatly appreciated.
Regards,
Luke Pretzie
  2 件のコメント
David Goodmanson
David Goodmanson 2016 年 10 月 27 日
Hi Luke, It's good to work out the details, but take a look at the 'conv' function too, info is in doc conv. The first case is just convolution of A with the vector [1 2 1]/4, and it's easy to generalize this to [1 2 3 2 1]/9 etc.
Image Analyst
Image Analyst 2016 年 10 月 27 日
Definitely agree. Use conv() instead of the loop you are doing.

サインインしてコメントする。

採用された回答

Star Strider
Star Strider 2016 年 10 月 27 日
You always need to take care in filter design so you get what you actually want. This demonstration code creates the FIR filters you described, and uses some Signal Processing Toolbox functions to analyse the filters and filter your signals:
[d,s,r] = xlsread('Luke Pretzie Copy of ECG data(2).xlsx');
t = d(:,1);
EKG = d(:,2);
N = 3;
v = ones(1,round(N/2));
b = conv(v, v);
a = sum(b);
figure(1)
freqz(b, a, 2^12)
EKG3 = filtfilt(b, a, EKG);
figure(2)
plot(t, EKG, t, EKG3,'-r')
grid
N = 5;
v = ones(1,round(N/2));
b = conv(v, v);
a = sum(b);
figure(3)
freqz(b, a, 2^12)
EKG5 = filtfilt(b, a, EKG);
figure(4)
plot(t, EKG, t, EKG5,'-r')
grid
If this does not give you the result you want, describe in some detail what you want. It’s probably possible to design a filter to do that.
  2 件のコメント
Luke Pretzie
Luke Pretzie 2016 年 10 月 31 日
Thank you very much, this filter works like a charm.
Regards, Luke Pretzie
Star Strider
Star Strider 2016 年 10 月 31 日
My pleasure.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeStatistics and Linear Algebra についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by