Contents

function [X, f, y, y2] = fftf(t, x, varargin)

% fftf - fft filter;
%   [X, f, y, y2] = fftf(t, x); with t the time vector and x the signal,
%   displays the original signal, the Fourier transform (absolute values)
%   and the reconstructed signal generated by the inverse transform ifft
%   with a selected subset of the frequencies.
%   By default, the frequencies in the filtered signal are cut at 1/8 the
%   sampling frequency.
%   The function returns X - reconstructed signal, f - vector of frequencies, y -
%   full vector of amplitudes, y2 - the filtered vector of amplitudes.
%   [X, f, y, y2] = fftf(t, x, cutoff); the user may set the cutoff
%   frequency in units of Hertz.
%   [X, f, y, y2] = fftf(t, x, cutoff, my_N); the user may select a number of my_N
%   amplitudes with highest absolute values to participate in the
%   reconstruction.
%   Examples
%   (suppose you have t and x defined already, both 1-dimensional vectors of the same length.)
%   fftf(t, x, 1e6); - reconstruct the signal with frequencies lower or
%   equal to cutoff value of 1MHz.
%   fftf(t, x, [], 20); - use only 20 biggest amplitudes for
%   reconstruction.
%   fftf(t, x, 1e6, 20); - select 20 biggest amplitudes within cutoff.
% Shmuel Ben-Ezra, Ultrashape ltd. August 2009

Verifying input

if ~any(size(t)==1),
   disp('Unexpected vector size! - should be 1D vectors.')
   return
end
if ~any(size(x)==1),
   disp('Unexpected vector size! - should be 1D vectors.')
   return
end
if length(t)~=length(x),
   disp('Unexpected vector size! - should be same length.')
   return
end

Definitions

Fs=1/(t(2)-t(1)); %sampling freq
N=length(x);
Nfft=2^nextpow2(N);
f=Fs/2*linspace(0,1,1+Nfft/2); % create freqs vector
cutoff_freq=Fs/8;
my_freqs=[];
if nargin>2,
    cutoff_freq=varargin{1};
end
if nargin>3,
    my_freqs=varargin{2};
end

main

y=fft(x,Nfft)/N; % perform fft transform
y2=filterfft(f, y, cutoff_freq, my_freqs); % filter amplitudes
%X=ifft(y2,'symmetric'); % the inverse transform. 'symmetric' is not recognized in older versions of matlab
X=ifft(y2); % inverse transform
X=X(1:N)*N;
ind1 = find(y2(1:1+Nfft/2)); % get the nonzero elements in y2
nf1 = length(ind1); % count nonzero elements

display

figname = 'fftf - FFT at work'; ifig = findobj('type', 'figure', 'name', figname); if isempty(ifig), ifig = figure('name', figname); % on my machine: ..., 'position', [360 120 600 800]); end figure(ifig); % first plot subplot(3,1,1) plot(t*1e6,x) xlabel('uSec') axis tight title('Original signal') % second plot subplot(3,1,2) yplot=abs(y(1:1+Nfft/2)); yplot=yplot/max(yplot); semilogy(f*1e-6, yplot, f(ind1)*1e-6, yplot(ind1), '.r'); xlabel('MHz') title('Amplitudes') legend('full spectrum', 'selected frequencies') % third plot subplot(3,1,3) plot(t*1e6,X) xlabel('uSec') if isempty(cutoff_freq), scutoff='No cutoff.'; else scutoff=sprintf('Cutoff = %g [Mhz]', cutoff_freq/1e6); end stitle3=sprintf('Reconstructed signal with %d selected frequencies; %s', nf1, scutoff); title(stitle3) axis tight return

function y2=filterfft(f, y, cutoff, wins)
nf=length(f);
ny=length(y);
if ~(ny/2+1 == nf),
    disp('unexpected dimensions of input vectors!')
    y2=-1;
    return
end

% cutoff filter
y2=zeros(1,ny);
if ~isempty(cutoff)
    ind1=find(f<=cutoff);
    y2(ind1) = y(ind1); % insert required elements
else
    y2=y;
end

% dominant freqs filter
if ~isempty(wins),
    temp=abs(y2(1:nf));
    y2=zeros(1,ny);
    for k=1:wins,   % number of freqs that I want
        [tmax, tmaxi]=max(temp);
        y2(tmaxi) = y(tmaxi); % insert required element
        temp(tmaxi)=0; % eliminate candidate from list
    end
end

% create a conjugate symmetric vector of amplitudes
for k=nf+1:ny,
    y2(k) = conj(y2(mod(ny-k+1,ny)+1)); % formula from the help of ifft
end
return

Published with MATLAB® 7.10