Problem: periodic noise in data (see blue curve below).
Sampling frequency: 40 Hz
Expected signal frequency: 10 Hz
Solution: low-pass filter (Butterworth)
Matlab already provides the framework for filter design (Butterworth filter in Matlab).
For explicit algorithms see http://www-users.cs.york.ac.uk/~fisher/mkfilter/
(specifically http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html)
Explanation for how filters work: http://www.electronics-tutorials.ws/filter/filter_8.html
Matlab filter usage: http://faculty.smu.edu/hgray/software/matlabfilterinstructions.pdf
Briefly:
Input is a 1D array of current measurements labeled "f". Output is array "filt".
Results:
1) No signal present. Expect a flat line in raw data.
Clearly measurement has periodic noise. Filter removes the noise very well. Note the ripple in the the beginning of the filtered data.
2) Signal is present around t = 1s.
Sampling frequency: 40 Hz
Expected signal frequency: 10 Hz
Solution: low-pass filter (Butterworth)
Matlab already provides the framework for filter design (Butterworth filter in Matlab).
For explicit algorithms see http://www-users.cs.york.ac.uk/~fisher/mkfilter/
(specifically http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html)
Explanation for how filters work: http://www.electronics-tutorials.ws/filter/filter_8.html
Matlab filter usage: http://faculty.smu.edu/hgray/software/matlabfilterinstructions.pdf
Briefly:
- the higher the order, the sharper the cutoff (see link above, 'how filters work')
- Nyquist frequency is equal to half of sampling frequency
- cannot apply filter to frequencies above Nyquist
Input is a 1D array of current measurements labeled "f". Output is array "filt".
f = f.*1E9; %convert to nA
time=(1:length(f)).*(1.0/40); % 40 Hz sampling rate
datahandle=plot(time,f);
xlabel('Time, s');
ylabel('Current, nA');
hold on;
sampling = 40; % sampling frequency = 40 Hz
nyquist = sampling/2; % nyquist frequency is half of the sampling frequency
cutoff = 15; % cutoff frequency = 15 Hz
n = 3; %filter order (higher is sharper cutoff)
Wn = cutoff/nyquist; % normalized cutoff frequency
[b,a]=butter(n,Wn); % get Butterworth coefficients
filt=filter(b,a,f); % apply filter to data
plot(time,filt,'Color','Red');
legend('Raw','Filtered');
Results:
1) No signal present. Expect a flat line in raw data.
Clearly measurement has periodic noise. Filter removes the noise very well. Note the ripple in the the beginning of the filtered data.
Zoom:
The spike in the beginning is an artifact of the algorithm. A quick solution is to insert fake data (repeated first value of the input array) in the beginning of the input array, apply the filter, and then remove the same amount of data from the beginning of the filtered array.
Matlab code:
f = f.*1E9; %convert to nA
time=(1:length(f)).*(1.0/40); % 40 Hz sampling rate
datahandle=plot(time,f,'Linewidth',2);
xlabel('Time, s');
ylabel('Current, nA');
hold on;
sampling = 40; % sampling frequency = 40 Hz
nyquist = sampling/2; % nyquist frequency is half of the sampling frequency
cutoff = 15; % cutoff frequency = 15 Hz
n = 3; %filter order (higher is sharper cutoff)
Wn = cutoff/nyquist; % normalized cutoff frequency
[b,a]=butter(n,Wn); % get Butterworth coefficients
pad_length = round(length(f)*0.1); % 10% of original data length
f = [ones(pad_length,1).*f(1); f]; %insert fake data in the beginning of the array
filt=filter(b,a,f); % apply filter to data
filt=filt(pad_length+1:length(filt)); %remove the fake data
plot(time,filt,'Color','Red','Linewidth',2);
legend('Raw','Filtered');
Results:
1) No signal present. The spike in the beginning is reduced.
2) Signal present around t = 1s.
No comments:
Post a Comment