About Butterworth filter design ?

Please; anyone can help me to design a Buttterworth filter for the attached data in the daily bases? i want to keep freq component up to 1.15^-8 Hz (1000 days and below). i take:

fc=1.15*10^-8;   % cuttoff freq (1000 day=1/1000*24*60*60)
fs=1.1*10^-5;   % sampling time (1 day)
% insert y data
[b,a]=butter(2,fc/(fs/2),'high' );
my = filter(b,a,y);plot(my)

Is that correct?

Answer 1

[Star Strider]

You have not clearly stated what you want to do.

You also have problems in your calculations, and your data are not uniformly sampled.

First, you need to use the Signal Processing Toolbox resample (link) function to regularly sample your data. Create a time vector for the resampling (interpolation) as:

tv = linspace(min(date), max(date), size(data,1));

From the documentation:

• y = resample(x,tx) resamples the values, x, of a signal sampled at the instants specified in vector tx. The function interpolates x linearly onto a vector of uniformly spaced instants with the same endpoints and number of samples as tx. NaNs are treated as missing data and are ignored.

Second, you are not calculating the sampling interval and therefore the sampling frequency correctly. Calculate them as:

Ts = mean(diff(tv));                                    % Sampling Interval (Days)
Fs = 1/Ts;                                              % Sampling Frequency (1/Day)
Fn = Fs/2;                                              % Nyquist Frequency (1/Day)

Then do a fft (link) on your resampled data to determine the frequencies you want to pass or reject, and specify your filter accordingly. I cannot do that because I do not know what you want to do.

Use the filtfilt function, not filter, to do the actual filtering.

[ahmad Saad]

thanks for your help Star Strider Actually, i want to remove trend and noise from the time series then i want to keep only the oscillations up to the fc=1.15*10^-8 Hz through a high pass filter , provided that the filtered data have the same length as the original data

[Star Strider]

My pleasure.

Other than the d-c offset (at 0 Hz), the lowest frequency you can identify in your data are about 1.97E-6 Hz. You can remove the d-c offset by subtracting the mean of your data, as I did here:

fidi = fopen('ahmad Saad data.txt', 'rt');
Read = textscan(fidi, '%s%f%f', 'HeaderLines',2, 'CollectOutput',1);
Time = Read{2}(:,1);
Data = Read{2}(:,2);
L = size(Time,1);
tv = linspace(min(Time), max(Time), L);
DataRsmp = resample(Data, tv);
Ts = mean(diff(tv));                                    % Sampling Interval (Days)
Fs = 1/Ts;                                              % Sampling Frequency (1/Day)
Fn = Fs/2;                                              % Nyquist Frequency (1/Day)

figure(1)
plot(tv, DataRsmp)
grid

DataDC0 = Data-mean(Data);
FTData = fft(DataDC0)/L;                                % Fourier Transform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                     % Frequency Vector (1/Day)
FvHz = Fv/(24*60*60);                                   % Frequency Vector (Hz)
Iv = 1:length(Fv);                                      % Index Vector

figure(2)
plot(Fv, abs(FTData(Iv))*2)
grid
axis([0  100    ylim])
xlabel('Frequency (1/Day)')
ylabel('Amplitude')

figure(3)
plot(FvHz, abs(FTData(Iv))*2)
grid
axis([0  0.0015    ylim])
xlabel('Frequency (Hz)')
ylabel('Amplitude')

In this context, ‘i want to keep only the oscillations up to the fc=1.15*10^-8 Hz through a high pass filter’ is a bit difficult for me to understand. Subtracting the mean may do everything you want.

The filtered signal using filtfilt (or filter) will be the same length as the original signal.