Strange result with Fourier transform after frequency removal

4 views (last 30 days)
I have a line profile that Im trying to obtain a "width measurement" on. The problem is that its polluted with wiggles.
I was thinking taking the FT, removing the highest component and then inverse FT would be a useful method.
However, my code is giving strange results. My x and y data is also attached
%look at FT to try and remove periodic ringing
F = fftshift(fft2(ydata)); %Take the FT
sF = log(abs(F)); %frequency spectrum
mxsF=max(sF(:)) %Determine max contributor
ind=find(sF==mxsF) %Find its frequency (i.e.. x-axis)
%View frequency current spectrum
subplot(1,4,4)
plot(sF,'b.')
grid on
hold on
%suppress peak freq
sF(ind)=0;
new=abs(fftshift(ifft(sF))); %inverse transform
plot(sF,'r-')
hold off
figure
plot(x,new)

Accepted Answer

Star Strider
Star Strider on 26 Jan 2016
The easiest way is to lowpass-filter your signal:
[d,s,r] = xlsread('Jason temp2.xlsx');
t = d(:,1);
s = d(:,2);
L = length(s);
Ts = mean(diff(t)); % Sampling Interval
Fs = 1/Ts; % Sampling Frequency
Fn = Fs/2; % Nyquist Frequency
FTs = fft(s)*2/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn; % Frequency Vector (Hz)
Iv = 1:length(Fv); % Index Vector
figure(1)
plot(Fv, abs(FTs(Iv)))
grid
Wp = 0.05/Fn; % Design Lowpass Filter
Ws = 0.15/Fn;
Rp = 1;
Rs = 30;
[n,Wn] = buttord(Wp,Ws,Rp,Rs);
[b,a] = butter(n,Wn);
[sos,g] = tf2sos(b,a);
figure(2)
freqz(sos) % Filter Bode Plot
s_f = filtfilt(sos, g, s); % Filter Signal
figure(3)
plot(t, s_f)
grid
  4 Comments
Jason
Jason on 26 Jan 2016
ahh. I don't have the buttord function as no signal processing toolbox. Is there a work around? thanks
Star Strider
Star Strider on 26 Jan 2016
My pleasure.
The best I can do in the absence of the Signal Processing Toolbox is Butterworth / Bessel / Chebyshev Filters from the University of York (U.K.).
The filter I used here has these coefficients:
b = [1.220454380697567e-03 4.881817522790266e-03 7.322726284185399e-03 4.881817522790266e-03 1.220454380697567e-03];
a = [1.000000000000000e+00 -2.897489668063421e+00 3.260788665023785e+00 -1.671326187500433e+00 3.275544606312290e-01];
Use those (and any you create from the York page) with the filter function. It’s not phase-neutral as is filtfilt, but it will likely work well enough.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!