I want to apply moving average filter! could you tell me how can I do that?

19 次查看(过去 30 天)
p = textread('Rainfall of Busan (11.5.1981 11.5.2021).txt');
dt = 1;
totalL = 14611;
T = [dt:dt:dt*totalL];
freqs = linspace(0,1/dt,totalL+1);
freqs = freqs(1:end-1);
y=fft(p);
I want to apply moving average filter! then what do I have to do...?
  1 个评论
Adam Danz
Adam Danz 2021-11-15
编辑:Adam Danz 2021-11-16
That's not specific enough. What kind of filter? How large is the window? Should the window move in increments of 1 or by window-width?

请先登录,再进行评论。

回答(2 个)

William Rose
William Rose 2021-11-16
@James James, everything @Adam Danz said is correct.
But here's an example, assuming you want a "flat" moving average:
t=0:100;
x=sin(2*pi*t/100)+rand(size(t))-.5; %for demonstration purposes
xf=zeros(size(x));
M=9; %window width, assumed to be odd
for i=(M+1)/2:101-(M-1)/2
xf(i)=sum(x(i-(M-1)/2:i+(M-1)/2))/M;
end
plot(t,x,'-r.',t,xf,'-b.',t,sin(2*pi*t/100),'-k');
legend('Raw','Filtered','No Noise');
Modify if you do not want output at every point. Modify as desired to handle the edges in a nicer way (e.g. extend the original signal to the sides before smoothing, with x(1) and x(end) respectively, or with a time-reversed and inverted extension, or whatever...). If you want a non-flat window, you have to do a bit more. For example, a triangular window of width M:
b=conv(ones(1,(M+1)/2),ones(1,(M+1)/2));
b=b/sum(b); %normalize so gain=1 at DC
Multipy elements of b by appropriate elements of x to get filtered version.

William Rose
William Rose 2021-11-16
Here's a method that hadles the edges better, and is more elegant, because it eliminates for loops.
t=0:100;
x=sin(2*pi*t/100)+rand(size(t))-.5; %for demonstration purposes
M=9; %window width, assumed to be odd
bf=ones(1,M)/M; %flat window
bt=conv(ones(1,(M+1)/2),ones(1,(M+1)/2))/((M+1)^2/4); %triangular window
xf=conv(x,bf);
xf=xf((M+1)/2:end-(M-1)/2);
xt=conv(x,bt);
xt=xt((M+1)/2:end-(M-1)/2);
plot(t,x,'-r.',t,xt,'-g.',t,xf,'-b.',t,sin(2*pi*t/100),'-k');
legend('Raw','Tri','Flat','No Noise');
Try.

产品

Community Treasure Hunt

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

Start Hunting!

Translated by