How to filter noise from time-frequency data and find natural frequency of a cantilever?

32 次查看（过去 30 天）
Rajesh
Rajesh 2016-5-11
评论： Image Analyst 2023-2-4
Hi,
I am studying the natural frequency of a cantilever. For my experiments I collect the deflection data (X and Y coordinates) of a mirror at the tip of the cantilever using a photodetector. I am using Fast Fourier Transform to convert the time domain data into frequency domain. I have attached two sample datasets. The code that I am using for the analysis is:
data=xlsread('sample_data_1.xlsx');
t=data(:,1);
y=data(:,3);
Fs=500;
nfft=2000;
Y=fft(y,nfft);
Y=Y(1:nfft/2);
mx=abs(Y);
f=(0:nfft/2-1)*Fs/nfft;
loglog(f,mx);
Ideally I should be getting a single peak which corresponds to the resonance/natural frequency of the cantilever. Surprisingly, I am getting multiple peaks with the peaks located at different frequencies for different dataset (some peaks coincide). This is due to noise from other sources. How can I get rid of the undesired frequencies? Any algorithm to filter out the noises? Is there an alternate way to analyse the dataset to obtain the natural frequency?
Regards,
Rajesh
  3 个评论
显示 1更早的评论
Rajesh
Rajesh 2016-5-11
编辑：Rajesh 2016-5-11
Currently I am observing the fluctuation in the cantilever position with no external force in effect. In the absence of an externally driving force the cantilever should vibrate at its own natural frequency. I am expecting noise from other sources because of mechanical and structural vibrations (other parts of the set up, vibration from the building etc.) present while I am doing the experiment. The cantilever I am studying (micro capillary) is quite sensitive to noise.
By variety of frequencies do you mean different modes of resonance? Yes that is possible but the frequencies should be integer multiples of each other. I don't see such a correlation.
I am observing multiple peaks with amplitude of oscillations in the range of 100s of microns. How do I figure out what is relevant and what isn't?
Image Analyst
Image Analyst 2023-2-4
Hmmm... your middle paragraph tells me that you don't have a good concept of Fourier (spectrum) analysis, and don't really know what the Fourier Transform means or how to interpret it. That's fine though, not everyone had a full semester course in it, like I did in grad school. But it would be good for you to at least learn the basics. I suggest you Google around and look for Fourier tutorials, especially those that talk about periodic signals.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Star Strider
Star Strider 2016-5-11
This is not a trivial problem. I took this approach and got the result in the plot. You will have to experiment to get the result you want, likely adding one or more additional sine terms to get a better fit. You will need the Signal Processing Toolbox to design and implement the filter. I used the core MATLAB function fminsearch. My objective function will work with the other curve fitting functions in the Statistics Toolbox and Optimization Toolbox. The fundamental natural frequency is 1/b(3) Hz. This should get you started, although you will need to experiment to get the result you want.
The code:
data=xlsread('Rajesh sample_data_2.xlsx');
t=data(:,1);
y=data(:,3);
Fs=500;
Fn = Fs/2; % Nyquist Frequency
L = size(t,1);
y = detrend(y); % Remove Linear Trend
fty = fft(y)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn; % Frequency Vector
Iv = 1:length(Fv); % Index Vector
figure(1)
plot(Fv, abs(fty(Iv)))
grid
Wp = 5/Fn; % Passband (Normalised)
Ws = 20/Fn; % Stopband (Normalised)
Rp = 1; % Passband Ripple
Rs = 25; % Stopband Ripple
[n,Wn] = buttord(Wp,Ws,Rp,Rs); % Filter Order
[b,a] = butter(n,Wn); % Transfer Function Coefficients
[sos,g] = tf2sos(b,a); % Second-Order-Section For Stability
figure(2)
freqz(sos, 2048, Fs)
yf = filtfilt(sos,g,y);
yu = max(yf);
yl = min(yf);
yr = (yu-yl); % Range of ‘y’
yz = yf-yu+(yr/2);
zt = t(yz(:) .* circshift(yz(:),[1 0]) <= 0); % Find zero-crossings
per = 2*mean(diff(zt)); % Estimate period
ym = mean(y); % Estimate offset
fit = @(b,x) b(1) .* exp(b(2).*x) .* (sin(2*pi*x./b(3) + 2*pi/b(4))) + b(5); % Objective Function to fit
fcn = @(b) sum((fit(b,t) - yf).^2); % Least-Squares cost function
s = fminsearch(fcn, [yr; -1; per; -1; ym]) % Minimise Least-Squares
xp = linspace(min(t),max(t));
figure(3)
plot(t,y,'b')
hold on
plot(t,yf,'y', 'LineWidth',2)
plot(xp,fit(s,xp), 'r', 'LineWidth',1.5)
hold off
grid
title('Sample Data 2')
xlabel('Time')
ylabel('Amplitude')
legend('Original Data', 'Lowpass-Filtered Data', 'Fitted Curve')
  2 个评论
uzzi
uzzi 2023-2-3
Hello,
I was also trying this code. It is working perfectly well until the line
s = fminsearch(fcn, [yr; -1; per; -1; ym])
where the error message is showing as
I also tried as
s = fminsearch(fcn, [double(yr); -1; per; -1; double(ym)])
but the error is still there.
Can you help me with this, please, @Star Strider?

请先登录，再进行评论。

Image Analyst
Image Analyst 2016-5-11
Depending on how it's moving, You may or may not see spikes. From what you say, I see no reason to expect it have spikes unless the structure is vibrating. If it's a sinusoid (in "resonance"), you'll see two spikes. If there are "harmonics" you will see multiple spikes as spacings a multiple of the first harmonic. But in general you will probably see a mountain shaped spectrum (due to the non-sinusoidal, wild and crazy motions that are present) with perhaps some small peaks on the sides of the mountain but only if it's vibrating. In other words, it is what it is. You got what you got. If you don't see an array of spikes, then that's just the way it is and you'll have to deal with it.

类别

在 Help Center 和 File Exchange 中查找有关 Statistics and Linear Algebra 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by