Help us understand the problem. What is going on with this article?

Matlabで時系列データをフィルタリングする

はじめに

周りの人に、この時系列データをいい感じにフィルタリングしたいんだけどMatlabでいい方法ない?ってきかれることありませんか?
たいていそういう人はバタワース関数を渡すと納得してくれるので、それをカットオフ周波数とサンプリング周波数を与えるだけに簡素化しました。

動作環境

Matlab R2018a
Signal Processing Toolboxが必要です。コマンドラインでverと打ってなければ使えません。

Ver .jpg

実行コード

clear all;
close all;
clc ;

%Data
Sampling = 1000; %Hz
t(:,1) = 0:1/Sampling:1;
Hz1 = 30;   %Hz
Hz2 = 150;  %Hz

ybase = 20*sin(Hz1*2*pi*t); %signal 
y1 = ybase + 1*sin(Hz2*2*pi*t)+ 5*rand(length(t),1) ; %Add Noise

% Designed Function
FilteredData_low = easyFilter(y1,Sampling,90);
FilteredData_band = easyFilter(y1,Sampling,[5 90]);
FilteredData_high = easyFilter(y1,Sampling,90,'high');

%result
subplot 411
plot(t,[ybase y1]);xlim([0.3 0.4]);

subplot 412
plot(t,[ybase FilteredData_low]);xlim([0.3 0.4]);
legend('lowpass')

subplot 413
plot(t,[ybase FilteredData_band]);xlim([0.3 0.4]);
legend('bandpass')

subplot 414
plot(t,[ybase ,FilteredData_high]);xlim([0.3 0.4]);
legend('highpass')

実行結果

時系列データy1に対して、90HzのLowパスフィルタ、5-90Hzのバンドパスフィルタ、90Hzのハイパスフィルタがかかります。赤線がもとの信号ybaseです。

easyFilter.jpg

ソースコード

使い方は以下です。数字はHz単位で統一してあります。

lowpass : FilteredData = easyFilter(y,1000,30);
bandpass : FilteredData = easyFilter(y,1000,[30 80]);
highpass : FilteredData = easyFilter(y,1000,30,'high');

使用関数の関係から、ナイキストの定理(1000Hzサンプリングに 2000HzのLowパス設定など)を超えるような設定すると駄目です。

function [FilteredData] = easyFilter(timeDomainData,Sampling_Hz,CutOffFrequency_Hz,varargin)
%EASYFILTER 時系列データに周波数指定でフィルタをかける
%  lowpass  : FilteredData = easyFilter(y,1000,30);
%  bandpass : FilteredData = easyFilter(y,1000,[30  80]);
%  highpass : FilteredData = easyFilter(y,1000,30,'high');
% 
%   Signal Processing Toolbox が必要

% y2_o3 wrote this file.
% 
% As long as you retain this notice you can do whatever you want
% with this stuff. If we meet some day, and you think this stuff
% is worth it, you can buy me a ** sushi ** in return.

if nargin == 3 % for low or bandpass filter
    [b,a] = butter(2,CutOffFrequency_Hz/(Sampling_Hz/2) );
    FilteredData = filtfilt(b,a,timeDomainData);

elseif nargin == 4 % for high pass filter
    [b,a] = butter(2,CutOffFrequency_Hz/(Sampling_Hz/2),varargin{1});
    FilteredData = filtfilt(b,a,timeDomainData);
else
    disp('Arguments Error@easyFilter: Quit');
end

end


解説(周囲対策用)

本フィルタは2次のバタワースフィルタをかけてあります。カットが緩やかな代わりに、
通過域の信号減衰がないないのが特徴です。
またfiltfilt関数を用いることで、リアルタイムでフィルタ処理をすると発生する
群遅延特性(位相遅れ)が解消されています。結果の図でlowpass後の信号のピークがずれていないのはそのためです。
オフライン解析では非常に有用な関数なので、ぜひお使いください。

おわりに

sushiライセンスを貼り付けました。だれかsushiをおごってくれませんかね。

hogemaccho
なんちゃってプログラマです。
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした