strange result from fft command

Hi,

I have written a code to calculate a derivative like

d/dx(c*df/dx)

with fft command (c is a function of x).

For example, if I define

f = sin(5*pi*x/L) and c(x)=x

I get correct result but if I define

f = sin(6*pi*x/L) and c(x)=x

I get wrong result!

In fact, whenever the coefficient of pi*x/L is an odd number I get correct result and when the coefficient of pi*x/L is an even number, the result is wrong.

changing c(x) does not affect this strange behavior.

If I use another function like f=cos(n*pi*x/L) , I get correct result when n is even and I do not get result when n is odd.

And finally, for other function of f, I do not get result.

Does anyone know the reason?

My Matlab version is R_2015b.

My code:

%%Constatnt Parameter %%

L=1.75;

h=L/150;

x=0:h:L;

f=-75:1:75; % frequency

w=2*pi*f/L; % wave number

%%calculating d/dx(c(x)*df/dx) in fourier space %%

n=7; % use odd and even number for n , and see the result

f=sin(n*pi*x/L); % defining function f(x), also test for f = cos(n*pi*x/L) and other function

y2=fftshift(fft(f)); % going to fourier space

y3=y2.*(1i.*w); % taking derivative of f(x) in fourier space

y4=ifft(ifftshift(y3)); % coming back to real space, the result is df/dx

c = x/L; % defining c(x)

y6 = c.*y4; % calculating c(x)*df/dx

y7=fftshift(fft(y6)); % going to fourier space

y8=y7.*(1i.*w); % taking derivative of c(x)*df/dx in fourier space

y9=ifft(ifftshift(y8)); % coming back to real space, the result is d/dx(c(x)*df/dx)

figure(1)

plot(x,y9,'.-b') % plot result from fft command

axis([0 L -(n*pi/L)^2 (n*pi/L)^2])

hold on

plot(x(2:150),diff(c(1:150).*diff(f)/h)/h,'r') % plot result from Finite difference Method

legend({'fft result','FD result'},'Location','northwest','FontSize',16)

Answer 1

[John BG]

Besides the point that you are applying FFT to an unstable signal, by this meaning, that it has no finite energy, cutting the signal and turning it periodic would help, there are different errors in your script.

1.- the way to apply fftshift and iffstshift. The following makes more sense, regarding the usage of fftshift and iffshift:

Y2=fftshift(f)
Y3=Y2.*(1i.*w)
y4=ifftshift(Y3)
c = x/L
y6 = c.*y4
Y7=fftshift(y6)
Y8=Y7.*(1i.*w)
y9=ifftshift(Y8)
figure(80);plot(x,y9,'.-b') ;grid on
axis auto

It reminds of the KEYSIGHT logo, doesn't it?

Since your time signal grows on both sides of the x axis, cannot have finite energy, therefore this cannot be the correct result, do you agree?

2.- fftshift and ifftshift useful to visualize, but you don't really need it here, just use fft ifft

Y20=fft(f)
Y3=Y2.*(1i.*w)
y4=ifft(Y3)
y6 = c.*y4
Y7=fft(y6)
Y8=Y7.*(1i.*w)
y9=ifft(Y8)
figure(9);plot(x,y9,'.-b') ;grid on

3.- Calculus derivative basics:

d/dx(c(x)*d/dx(f(x))) = cdot*fdot + c*fdotdot

where

cdot = d/dx(c(x))

fdot = d/dx(f(x))

fdotdot = d/dx(d/dx(f(x)))

in MATLAB:

cdot=ifft(1j.*w.*fft(c))
fdot=ifft(1j.*w.*fft(f))
fdotdot=-ifft(w.*fft(f))

So one way to start solving this problem would be:

cdot*fdot+c*fdotdot = IFFT(FFT(cdot*cdot+c*fdotdot)) = IFFT(FFT(cdot*cdot))+IFFT(FFT(c*fdotdot)))

Or perhaps you would like to start with a more compact expression:

cdot*fdot+c*fdotdot = IFFT(jw*FFT(c))*IFFT(jw*FFT(c))-c*IFFT(w*FFT(f)) =

-c*IFFT(w*FFT(f)) + j*IFFT(w*FFT(c))^2

3.- Perhaps you want to consider taking a shorter route: time domain

y20=diff(f)
c(end)=[]
c.*y20
y23=c.*y20
y24=diff(y23)
y24=[y24 y24(end)]
plot([1:1:length(y20)-2],y20(1:end-2),'r',[1:1:length(y23)-1],y23(1:end-1),'g',[1:1:length(y24)],y24,'b');
grid on
figure(20);plot(x(2:150),diff(c(1:150).*diff(f)/h)/h,'r');
grid on

.

If you find this answer useful, please click on thumbs-up vote above, thanks in advance

John

[e-mail]

Answer 2

[John BG]

Good, you paid attention to the last part, where I did same mistake I mentioned to avoid earlier on.

Let me correct time domain answer with the following:

% input signals
c = x/L
f=sin(n*pi*x/L)

% apparent solution: diff(c(1:150).*diff(f)/h)/h
y20=diff(f)
c(end)=[]
c.*y20
y23=c.*y20
y24=diff(y23)
y24=[y24 y24(end)]
plot([1:1:length(y20)-2],y20(1:end-2),'r',[1:1:length(y23)-1],y23(1:end-1),'g',[1:1:length(y24)],y24,'b');grid on
figure(20);plotyy([2:1:150],f(2:150),[2:1:150],diff(c(1:150).*diff(f)/h)/h);grid on;grid minor;axis auto

 applying d/dx(c(x)*d/dx(f(x)))=cdot*fdot+c*fdotdot

c2=c
f2=f
c2dot=diff(c2)/h
f2dot=diff(f2)/h

c2(end)=[]                 % making sure all vector to combine have same length
f2(end)=[]

f2dotdot=diff(f2dot)/h   % making sure all vector to combine have same length

c2(end)=[]
f2(end)=[]
c2dot(end)=[]
f2dot(end)=[]

y=c2dot.*f2dot+c2.*f2dotdot

nc2=[1:1:length(c2)]
% nf2, nc2dot, nf2dot, nf2dotdot are exactly same as nc2, 
% so using nc2 as reference for all signals

figure(1);plot(nc2,c2,nc2,f2,nc2,y)

figure(2);plotyy(nc2,f2,nc2,y)
grid on;grid minor

Now, if you have them off, turn on the following variable descriptors: Range, Min,Max, Mean, Var, Std in MATLAB workspace

y and y24 look alike, but y is the correct output signal.

Now

1.- Regarding you quoting some one that suggests in book loss of 'product rule' because of sampling, and not preserving Hermitian property of the operator (which one, the product?) .. could you please quote book title, ISBN and author(s) so we the readers of your question can read the source?

Is it possible for you to detail the area of application of this equation? is it civil engineering? trying to bring down something? or trying to make something robust enough not to be brought down? is is electronics?

2.- You have chosen a sampling step h, not the readers of your question. I can only work with the data you have supplied. If you want you can increase accuracy by reducing the step size, increasing window of observation of x, f, .. but you, not the readers of your question, chose to work with samples, and then you mention loss of, product rule? of Hermitian operator? What is it you are after, the correct output of your system?

what matrix is it that you want to preserve along the processing of the input signal? if you don't show the bigger picture, it is difficult to help.

3.- we need ..? and click a web link? i don't need to use fft shift, you do. Or you think you do. And it would be kindly appreciated defining the question concisely without waiting to a second round to ask for clicks on literature references that you don't want to mention to the readers of your question. Please define the question correctly, quote book reference details.

4.- Let me show you the symbolic result with MuPAD:

d/dx(c(x))=1/L is a constant, w is a constant, c'f'+cf'' goes f'+cf'' do you agree?

if you really want to preserve Continuous Time properties, then don't use FFT, use DFT or DCT with the right parameters.

FFT is a particular case of DFT that trades off accuracy for processing rate.

Let be F=FT(f) C=FT(c) then FT(f'+cf'') = jwF + FT(cf'')

Do you still really think you must transform to frequency and then back to time to get the function that MuPAD concisely shows?

Then go ahead, write a script that performs

-c*IFT(w*FT(f)) + j*IFT(w*FT(c))^2

and, then, if you want, let us know if it works. The working script would also be appreciated.

If you think that Discrete Time domain is going to deteriorate something that you hold dearly, then don't use it, work in symbolic, but showing up with DT question and then arguing that sampling brings it all down, really, it's disappointing. Make a choice. To be or not to be .. Hamlet?

5.- Just in case not yet noticed, you chose exactly same variable name for the frequency reference vector f and the input function, in my opinion unfortunate choice. Try renaming the input function to something different

Good luck, really, keep us posted.

John [e-mail]