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

51 views (last 30 days)
Junu Lee on 2 Sep 2020
Commented: Star Strider on 4 Oct 2020
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.
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.

Star Strider on 2 Sep 2020
I am not certain what you want to do.
Try this:
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)
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.
##### 2 CommentsShowHide 1 older comment
Star Strider on 4 Oct 2020
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:
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)
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. 