How can I obtain the frequency response and the cutoff frequency from sensor data text file?

I use the laser sensor for measuring the distance to obtain data.

So I obtain a text file what I attach.

The first columb is the distance data measured by laser sensor.

The period of the sensor is 10ms.

I make a code time to data as below.

data = load('Suctiontest_4.txt');
A = data(:, 1);
AA = A';
t = 0:0.01:35.81
plot(t, AA)
grid on

Then I wanna obtain the cutoff frequency because there is noise in this data.

How can I get the cutoff frequency and how can I plot the bode?

And how can I get the transfer function of this sensor?

Thank you.

Answer 1

[Star Strider]

I am not certain what you want to do.

Try this:

data = load('Suctiontest_4.txt');
AA = data(:, 1);
% AA = A';
t = 0:0.01:35.81;
plot(t, AA)
grid

Fs =  100;
Fn = Fs/2;
L = size(data,1);
datam = mean(data);
FTdata = fft(data - datam)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTdata(Iv,1))*2)
grid
xlim([0 5])

Rp =  1;
Rs = 60;
Wp = 0.5/Fn;
Ws = 0.6/Fn;
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);
[z,p,k] = ellip(n,Rp,Rs,Wp,'low');
[sos,g] = zp2sos(z,p,k);

figure
freqz(sos,2^14,Fs)

AA_filt = filtfilt(sos, g, AA);

figure
plot(t, AA, '-b')
hold on
plot(t, AA_filt, '-r', 'LineWidth',1.5)
hold off
grid

AAm = mean(AA_filt);
FTAA = fft(AA_filt - AAm)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:numel(Fv);

figure
subplot(2,1,1)
plot(Fv, 20*log10(abs(FTAA(Iv))))
grid
xlim([0 5])
ylabel('dB')
subplot(2,1,2)
plot(Fv, rad2deg(unwrap(angle(FTAA(Iv)))))
grid
xlim([0 5])
xlabel('Frequency (Hz)')
ylabel('Phase (°)')

Experiment with the filter passband (‘Wp’) and stopband (‘Ws’) frequencies (here 0.5 and 0.6 Hz respectively) to get the result you want. Note that for any lowpass filter, Ws > Wp.

[Junu Lee]

Could you tell me how i obtain the cut off frequency by -3dB?

[Star Strider]

The filter cutoff frequency ‘Wp’ is whatever you define it to be. Here, I used the Fourier transform of the signal (figure(2) in my code) to determine that, and defined it as 0.5 Hz. It and the stopband freqeuncy ‘Ws’, can be whatever you want them to be, between 0 and the Nyquist frequency, excluding both of those limits, so the normalised frequencies are between (but not including) 0 and 1. For a lowpass filter such as this one, the stopband frequency ‘Ws’ must be greater than the passband frequency, ‘Wp’. Here, I set ‘Ws’ to 0.6 Hz to give a relatively short transition region. That still results in a short, efficient filter.

The corresponding -3 dB and -6 dB points from the data is a bit more complicated.

Add these to my previous coed:

dB36 = 10.^([-3 -6]/20);                                                            % Define -3dB & -6dB Values
[FTmax,ix] = max(abs(FTdata(Iv,:))*2);                                              % Maximum For Each Column
for k = 1:size(FTdata,2)
    FTv = abs(FTdata(Iv,k))*2;                                                      % Fourier Transformm Of Each Column
    minv = find(islocalmin(FTv, 1));                                                % Find Local Minima
    min1(k) = minv(1);                                                              % Define First Minimum
    Fv1dB36(k,:) = interp1(FTv(ix(k)-1:ix(k)),Fv(ix(k)-1:ix(k)),dB36*FTmax(k));     % Interploate To Find -3dB & -6dB Points
    Fv2dB36(k,:) = interp1(FTv(ix(k):min1(k)),Fv(ix(k):min1(k)),dB36*FTmax(k));     % Interploate To Find -3dB & -6dB Points
end

figure
plot(Fv, abs(FTdata(Iv,1))*2)
hold on
plot(Fv2dB36(1,:), dB36*FTmax(1), 'r+')
hold off
grid
xlim([0 5])
text(Fv2dB36(1,:), dB36*FTmax(1), compose('\\leftarrow %2d dB %.3f Hz', [-3 -6; Fv2dB36(1,:)]'), 'HorizontalAlignment','left', 'VerticalAlignment','middle')

so the full code is now:

data = load('Suctiontest_4.txt');
AA = data(:, 1);
% AA = A';
t = 0:0.01:35.81;

figure
plot(t, AA)
grid
xlabel('t')
ylabel('AA')

Fs =  100;
Fn = Fs/2;
L = size(data,1);
datam = mean(data);
FTdata = fft(data - datam)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:numel(Fv);

dB36 = 10.^([-3 -6]/20);                                                            % Define -3dB & -6dB Values
[FTmax,ix] = max(abs(FTdata(Iv,:))*2);                                              % Maximum For Each Column
for k = 1:size(FTdata,2)
    FTv = abs(FTdata(Iv,k))*2;                                                      % Fourier Transformm Of Each Column
    minv = find(islocalmin(FTv, 1));                                                % Find Local Minima
    min1(k) = minv(1);                                                              % Define First Minimum
    Fv1dB36(k,:) = interp1(FTv(ix(k)-1:ix(k)),Fv(ix(k)-1:ix(k)),dB36*FTmax(k));     % Interploate To Find -3dB & -6dB Points
    Fv2dB36(k,:) = interp1(FTv(ix(k):min1(k)),Fv(ix(k):min1(k)),dB36*FTmax(k));     % Interploate To Find -3dB & -6dB Points
end

figure
plot(Fv, abs(FTdata(Iv,1))*2)
hold on
plot(Fv2dB36(1,:), dB36*FTmax(1), 'r+')
hold off
grid
xlim([0 5])
text(Fv2dB36(1,:), dB36*FTmax(1), compose('\\leftarrow %2d dB %.3f Hz', [-3 -6; Fv2dB36(1,:)]'), 'HorizontalAlignment','left', 'VerticalAlignment','middle')

Rp =  1;
Rs = 60;
Wp = 0.5/Fn;
Ws = 0.6/Fn;
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);
[z,p,k] = ellip(n,Rp,Rs,Wp,'low');
[sos,g] = zp2sos(z,p,k);

figure
freqz(sos,2^14,Fs)

AA_filt = filtfilt(sos, g, AA);

figure
plot(t, AA, '-b')
hold on
plot(t, AA_filt, '-r', 'LineWidth',1.5)
hold off
grid

AAm = mean(AA_filt);
FTAA = fft(AA_filt - AAm)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:numel(Fv);

figure
subplot(2,1,1)
plot(Fv, 20*log10(abs(FTAA(Iv))))
grid
xlim([0 5])
ylabel('dB')
subplot(2,1,2)
plot(Fv, rad2deg(unwrap(angle(FTAA(Iv)))))
grid
xlim([0 5])
xlabel('Frequency (Hz)')
ylabel('Phase (°)')

This displays one signal. Putting the figure(2) call in a for loop and looping over the columns would display all of them. The corresponding -6 dB and -3 dB levels on the rising part of the peak are the same for all the columns: 0.0198 and 0.0140 Hz respectively, so I did not include them in the plot.