Extract time signal of specific frequences from a given time series signal

Hi all,

I want to extract time signal of two specific frequencies and their combinations from a given signal.

I tried the following codes and I found the extracted signals are very differnent from the corresponding ones.

Is this filter design correct? If not, can you tell how to extract the correct signals? Thank you.

figure;
fs = 1000; % sampling frequency 1000 [Hz]
T = .1; % signal total Time 0.1 [s]
dt = 1/fs; % time resolution
t = dt:dt:T; % time axis

speed = 80; % rotating speed

x1 = 10*sin(2*pi*speed*t+10);
x2 = 3*sin(2*pi*(2*speed)*t+10);
noise = 0.05 * rand(1,numel(t));
x = x1 + x2 + noise;


subplot(3,2,1);
plot(t,x1);
title('original x1');

subplot(3,2,3);
plot(t,x2);
title('original x2');

subplot(3,2,5);
plot(t,x);
title('original x');

% filtered x1
[z,p,k] = butter(4,[speed-1, speed+1]/(fs/2),'bandpass');
[sos,g] = zp2sos(z,p,k); %sos representation of the filter
filteredX1 = filtfilt(sos,g,x);
subplot(3,2,2);
plot(t,filteredX1);
title('filtered x1');


% filtered x2
[z,p,k] = butter(4,[2*speed-1, 2*speed+1]/(fs/2),'bandpass');
[sos,g] = zp2sos(z,p,k); %sos representation of the filter
filteredX2 = filtfilt(sos,g,x);
subplot(3,2,4);
plot(t,filteredX2);
title('filtered x2');

% combined x1 + x2
combined = filteredX1 + filteredX2;
subplot(3,2,6);
plot(t,combined);
title('combined of filtered x1 + x2');  

Answer 1

[Paul]

Hi Chenghao,

I suspect the issue is related to the duration of x relative to the dynamics of the transient response of the filter.

fs = 1000; % sampling frequency 1000 [Hz]

Extend the duration of the data

%T = .1; % signal total Time 0.1 [s]
T = 5;
dt = 1/fs; % time resolution
t = dt:dt:T; % time axis
speed = 80; % rotating speed
x1 = 10*sin(2*pi*speed*t+10);
x2 = 3*sin(2*pi*(2*speed)*t+10);
noise = 0.05 * rand(1,numel(t));
x = x1 + x2 + noise;

[z,p,k] = butter(4,[speed-1, speed+1]/(fs/2),'bandpass');
[sos,g] = zp2sos(z,p,k); %sos representation of the filter
[b,a]   = zp2tf(z,p,k);

Use filter, instead of filtfilt to illustrate

%filteredX1 = filtfilt(sos,g,x);
filteredX1 = filter(b,a,x);
figure
plot(t,filteredX1)

It takes about 2 seconds for the filtered output to reach steady state. filtfilt does attempt to minimize transients, but I'm not sure what to expect if the duration of the input is so much shorter than the "duration" of the filter's transient response.