Differentiation from FFT finding extrema

96 views Asked by At

I'm trying to find zeros of a function. See my code below.

Because fft expects a numerical array, I didn't define the symbolic function to use fzero.

However, this approach is not accurate and depend on step. Do you have a better idea?

step=2000;
x=0:pi/step:2*pi;
y= 4+5*cos(10*x)+20*cos(40*x)+cos(100*x);
fy = fft(y');
fy(1:8) =0;
fy(12:end-10)=0; 
fy(end-6:end)=0;
ffy = ifft(fy);
t=diff(ffy);
x=0:pi/step:2*pi-pi/step;
plot(x,t)
indices= find(t<5e-4 & t>-5e-4);
1

There are 1 answers

3
eigenchris On BEST ANSWER

You could proceed along the array t and look for points where the values change sign. That would indicate the presence of a zero.

Actaully, MATLAB's fzero function uses a similar method. You said you didn't use it because you required an array, rather than an anonymous function, but you could convert the array into an anonymous function using simple linear interpolation like so:

func = @(k) interp1(x,t,k);  % value of t(x) interpolated at x=k
fzero(func,initial_value);


EDIT : Just to clarify what I mean. If you have an array t and you want to find its zeros...

f = 5;                            % frequency of wave in Hz
x = 0:0.01:1;                     % time index
t = cos( 2*pi*f*x );              % cosine wave of frequency f

zeroList = [];                    % start with an empty list of zeros
for i = 2:length(t)               % loop through the array t

    current_t = t(i);             % current value of t
    previous_t = t(i-1);          % previous value of t

    if current_t == 0                  % the case where the zero is exact
        newZero = x(i);
        zeroList = [zeroList,newZero];
    elseif current_t*previous_t < 0    % a and b have opposite sign if a*b is -ve
        % do a linear interpolation to find the zero (solve y=mx+b)
        slope = (current_t-previous_t)/(x(i)-x(i-1));
        newZero = x(i) - current_t/slope;
        zeroList = [zeroList,newZero];
    end

end

figure(1); hold on;
axis([ min(x) max(x) -(max(abs(t))) (max(abs(t))) ]); 
plot(x,t,'b');
plot(x,zeros(length(x)),'k-.');
scatter(zeroList,zeros(size(zeroList)),'ro');

The zeros I get are correct:

enter image description here