Hi all, 
  I am using a frequency operator to act as a differentiator like

  i*w <==> d/dt
here w is the frequency defined in fourier space.

I have my code written in matlab as follow

n   = 4;        % dx = 1/2^n;
Ld  = 100;      % the bound value of the range
TNT = Ld*2^n;
NX  = 2.^(ceil(log(TNT)/log(2))); % number of sample
dx  = Ld/NX;;                     % sampling interval
x   = (-Ld/2:dx:Ld/2-dx);         
tnn = 0:(NX/2); 
nn  = [tnn  -fliplr(tnn(2:end-1))]; 
w   = 2*pi*nn/(NX*dx);
% I take the frequency [0 1 2 ... (NX/2) -(NX/2-1) -(NX/2-2) ... -1]

y  = sin(x);
dy = ifft(i*w.*fft(y));
subplot(2,1,1);
plot(x, y);
axis([-10 10 0 1]);
grid on;
subplot(2,1,2);
plot(x, dy);
axis([-10 10 0 1]);
grid on;

The resulting plot shows that there is a slight phase shift. I guess
it's due to the error of fft(i.e the non-zero imaginary after taking
the inverse-fft, is that right?) By plotting the amplitude of the
result, the phase shift vanishes

subplot(2,1,1);
plot(x, abs(y));
axis([-10 10 0 1]);
grid on;
subplot(2,1,2);
plot(x, abs(dy();
axis([-10 10 0 1]);
grid on;

By the way, I wonder why there is a sharp bounce at the boundary which
introduces an error. In my code, it seems that Ld is a significant
parameter. Decreasing Ld will lead to the instability of dy (some
oscillating wave will be added to the result). What's the reason? It
is so weird. I think the sampling interval will only affected by the
value of n. How can I choose n and Ld to balance the efficiency and
accuracy?

Thanks in advance.

M...@21cn.com (Hakuna M. C.) wrote in message news:<3...@posting.google.com>...
> Hi all, 
>   I am using a frequency operator to act as a differentiator like
> 
>   i*w <==> d/dt
> here w is the frequency defined in fourier space.

Almost. What you state is valid for continuous functions. With matlab
(and any other FFT routine) you work with discrete sequences. You need 
to use the differentiation theorem of the z transform. Very crudely, 
you get something like 

   "dx(n-1)/dt" = 2*real(IFFT(k*X(k))), k=0,...,N/2

where dx(n)/dt is a "quasi differential" of the continuous x(t).
Note the circular time shift!    

DISCLAIMER: I have just found out about these things, as I'm working 
with a problem that involves carrying "quasi differentials" of discrete 
numbers through the DFT/IDFT. There may be (and probably is) a few 
wrong details here, but the difference between the Fourier transform 
and z transform is cruical. I'd appreciate any input on the Hilbert 
transform (which is not quite the same as the Hilbert transformer).

Rune

"Rune Allnor" <a...@tele.ntnu.no> wrote in message
news:f...@posting.google.com...
> M...@21cn.com (Hakuna M. C.) wrote in message
news:<3...@posting.google.com>...
> > Hi all,
> >   I am using a frequency operator to act as a differentiator like
> >
> >   i*w <==> d/dt
> > here w is the frequency defined in fourier space.
>
> Almost. What you state is valid for continuous functions. With matlab
> (and any other FFT routine) you work with discrete sequences. You need
> to use the differentiation theorem of the z transform. Very crudely,
> you get something like
>
>    "dx(n-1)/dt" = 2*real(IFFT(k*X(k))), k=0,...,N/2
>
> where dx(n)/dt is a "quasi differential" of the continuous x(t).
> Note the circular time shift!
>
> DISCLAIMER: I have just found out about these things, as I'm working
> with a problem that involves carrying "quasi differentials" of discrete
> numbers through the DFT/IDFT. There may be (and probably is) a few
> wrong details here, but the difference between the Fourier transform
> and z transform is cruical. I'd appreciate any input on the Hilbert
> transform (which is not quite the same as the Hilbert transformer).

Rune,

In O&S 1st ed. Chapter 7 they say that the "Hilbert Transform Relations" are
(I paraphrase greatly here) simply things regarding observations about
causality, minimum phase (or not), real parts, etc. etc.  So, I wonder if
that's not all there is to it?

Fred

Fred

"Hakuna M. C." <M...@21cn.com> wrote in message
news:3...@posting.google.com...
> Hi all,
>   I am using a frequency operator to act as a differentiator like
>
>   i*w <==> d/dt
> here w is the frequency defined in fourier space.
>
> I have my code written in matlab as follow
>
> n   = 4;        % dx = 1/2^n;
> Ld  = 100;      % the bound value of the range
> TNT = Ld*2^n;
> NX  = 2.^(ceil(log(TNT)/log(2))); % number of sample
> dx  = Ld/NX;;                     % sampling interval
> x   = (-Ld/2:dx:Ld/2-dx);
> tnn = 0:(NX/2);
> nn  = [tnn  -fliplr(tnn(2:end-1))];
> w   = 2*pi*nn/(NX*dx);
> % I take the frequency [0 1 2 ... (NX/2) -(NX/2-1) -(NX/2-2) ... -1]
>
> y  = sin(x);
> dy = ifft(i*w.*fft(y));
> subplot(2,1,1);
> plot(x, y);
> axis([-10 10 0 1]);
> grid on;
> subplot(2,1,2);
> plot(x, dy);
> axis([-10 10 0 1]);
> grid on;
>
> The resulting plot shows that there is a slight phase shift. I guess
> it's due to the error of fft(i.e the non-zero imaginary after taking
> the inverse-fft, is that right?) By plotting the amplitude of the
> result, the phase shift vanishes
>
> subplot(2,1,1);
> plot(x, abs(y));
> axis([-10 10 0 1]);
> grid on;
> subplot(2,1,2);
> plot(x, abs(dy();
> axis([-10 10 0 1]);
> grid on;
>
> By the way, I wonder why there is a sharp bounce at the boundary which
> introduces an error. In my code, it seems that Ld is a significant
> parameter. Decreasing Ld will lead to the instability of dy (some
> oscillating wave will be added to the result). What's the reason? It
> is so weird. I think the sampling interval will only affected by the
> value of n. How can I choose n and Ld to balance the efficiency and
> accuracy?

In effect, you're trying to apply a FIR filter.  Have you looked at the
differentiator FIR filters that can be designed using the Parks-McClellan
program?  Seems that might lend some insight - as perfection isn't
attainable.  There's a good treatment of the subject in Rabiner & Gold
"Theory and Application of Digital Signal Processing" p.164 section 3.36
where they discuss the difference between Type 3 and Type 4 differentiators.
A typical approach is to have the error grow with frequency.  Your approach
assumes that a perfect differentiator is possible.  You might calculate the
ifft of your differentiator frequency function and see what kind of impulse
response it has.

As far as your implementation goes, you're multiplying in frequency which
has a corresponding convolution operation in time.  It doesn't look like
you've taken temporal aliasing into account.  The rule is this:
in order to convolve two arrays in time / multiply in frequency / with the
starting arrays of length M and N, then the transformed length needs to be
at least M+N-1 with sufficient zeros to not have the leading and trailing
edges of the convolution overlap and alias.  You need to append at least N-1
zeros to the array of length M and at least M-1 zeros to the array of length
N - with the end result that both arrays are the same length of course.  Or,
if M and N are equal, you need to add at least N-1 zeros to both arrays.

For the kind of analysis you're doing, you may wish to append the zero
samples in such a way as to keep the center of the waveform at t=0.  This
means the array will have the zeros appended in the middle and will have
nonzero values at t=0 and above and at t=(N-1)*tau and below where (N-1)*tau
is the length of the temporal epoch that your N samples cover with tau being
the temporal sample interval.  So, you may want the zeros to be centered
around t=(N-1)*tau/2.   This keeps the signal centered at zero - so there's
no delay.  You can also do this with the differentiator impulse response.

Note that doing the fft assumes that the time function is periodic so the
first value is assumed to repeat at t=N*tau.

In the end:  you can either convolve in time or multiply in frequency.  It
is just more computationally efficient to fft/multiply/ifft.  Otherwise you
can use "conv".  The results in the coresponding time and frequency domains
should be exactly the same and usually are very close in the type of thing
your program represents.  So, don't be looking so quickly for numerical
errors.  You need to convolve in time with something that makes sense and
with enough zeros in the arrays so the results of convolution don't overlap.
If you use "conv", the results are extended for you.  But, if you use
fft/multiply/ifft or "circular convolution" then this isn't done for you -
you have to do it yourself.  That is you may correspondingly multiply in
frequency with something that makes sense as above.  i*w isn't one of those
I believe.

Fred

"Rune Allnor" <a...@tele.ntnu.no> wrote in message
news:f...@posting.google.com...
> M...@21cn.com (Hakuna M. C.) wrote in message
news:<3...@posting.google.com>...
> > Hi all,
> >   I am using a frequency operator to act as a differentiator like
> >
> >   i*w <==> d/dt
> > here w is the frequency defined in fourier space.
>
> Almost. What you state is valid for continuous functions. With matlab
> (and any other FFT routine) you work with discrete sequences. You need
> to use the differentiation theorem of the z transform. Very crudely,
> you get something like
>
>    "dx(n-1)/dt" = 2*real(IFFT(k*X(k))), k=0,...,N/2
>
> where dx(n)/dt is a "quasi differential" of the continuous x(t).
> Note the circular time shift!
>
> DISCLAIMER: I have just found out about these things, as I'm working
> with a problem that involves carrying "quasi differentials" of discrete
> numbers through the DFT/IDFT. There may be (and probably is) a few
> wrong details here, but the difference between the Fourier transform
> and z transform is cruical.

Rune,

What, would you say is the crucial difference?  My notion (without thinking
about it) is that the Discrete Fourier Transform computed to yield a
discrete spectral transform is evaluated on a unit circle and the Z
Transform need not be but usually is anyway??

Or is it that the z transform is analogous to the Laplace transform and
systems need to be causal??

Fred

"Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message
news:<u...@centurytel.net>...
> "Rune Allnor" <a...@tele.ntnu.no> wrote in message
> news:f...@posting.google.com...

> > I'd appreciate any input on the Hilbert
> > transform (which is not quite the same as the Hilbert transformer).
> 
> Rune,
> 
> In O&S 1st ed. Chapter 7 they say that the "Hilbert Transform Relations" are
> (I paraphrase greatly here) simply things regarding observations about
> causality, minimum phase (or not), real parts, etc. etc.  So, I wonder if
> that's not all there is to it?

Well, yes... but I'm a bit curious about the full impact of those 
observations/properties. If a real-valued sequence is causal, the real 
and imaginary parts of its spectrum must be Hilbert transform pairs. 
If you know, say, the real part of the spectrum you can compute the 
imaginary part. And vice versa. To me, this hints that the Hilbert 
tranform could be quite a powerful tool, if handeled correctly. 

Suppose there are other Hilbert transform pairs, like perhaps the 
magnitude and phase components of spectra. Suppose there are other 
variations of the Hilbert transform. The Fourier-Hankel tranform is, 
for instance, completely analogous to the 2D FT but in a cylindrical 
coordinate system. I am getting curious about what the Hilbert 
transform is, first of all, and then if it can be generalized in 
one way or another. 

If the Hilbert transform could be generalized to allow you to state 
constraints on general formulations of signals in order for the 
signals to obtain "arbitrary" properties, I certainly would like to 
know about it.

Rune

"Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message
news:<w...@centurytel.net>...
> "Rune Allnor" <a...@tele.ntnu.no> wrote in message
> news:f...@posting.google.com...
> > M...@21cn.com (Hakuna M. C.) wrote in message
>  news:<3...@posting.google.com>...
> > > Hi all,
> > >   I am using a frequency operator to act as a differentiator like
> > >
> > >   i*w <==> d/dt
> > > here w is the frequency defined in fourier space.
> >
> > Almost. What you state is valid for continuous functions. With matlab
> > (and any other FFT routine) you work with discrete sequences. You need
> > to use the differentiation theorem of the z transform. Very crudely,
> > you get something like
> >
> >    "dx(n-1)/dt" = 2*real(IFFT(k*X(k))), k=0,...,N/2
> >
> > where dx(n)/dt is a "quasi differential" of the continuous x(t).
> > Note the circular time shift!
> >
> > DISCLAIMER: I have just found out about these things, as I'm working
> > with a problem that involves carrying "quasi differentials" of discrete
> > numbers through the DFT/IDFT. There may be (and probably is) a few
> > wrong details here, but the difference between the Fourier transform
> > and z transform is cruical.
> 
> Rune,
> 
> What, would you say is the crucial difference?  My notion (without thinking
> about it) is that the Discrete Fourier Transform computed to yield a
> discrete spectral transform is evaluated on a unit circle and the Z
> Transform need not be but usually is anyway??
> 
> Or is it that the z transform is analogous to the Laplace transform and
> systems need to be causal??
> 
> Fred

I need some maths to explain the problem.

Define the z tranform (ZT) as

                        inf
     X(z) =  Z{x(n)} =  sum    x(n)z^(-n)                     [1]
                       n=-inf

and DFT as

                          inf
     X(w) =  DFT{x(n)} =  sum   x(n)exp(-jwn).                [2]
                         n=-inf

===========================================================================

Properties of the Tranforms:

Time Shift:

                  inf
    Z{x(n-m)} =   sum   x(n-m)z^(-n)                          [3.a]
                n=-inf

                   ( p=n-m =>  n=p+m )

                  inf
              =   sum    x(p)z^(-p-m)                         [3.b]
                 p=-inf

              =  z^(-m) Z{x(n)}                               [3.c]


----------------------------------------------------------------------

                   inf
    DFT{x(n-m)} =  sum     x(n-m) exp(-jwn)                   [4.a]
                  n=-inf
      
                   inf
                =  sum     x(p) exp(-jw(p+m))                 [4.b]
                  n=-inf
             
                = exp(-jwm) DFT {x(n)}.                       [4.c]


=========================================================================

Differentiation if w/z domain:

    dZ{x(n)}        d   inf
    --------    = ---   sum  x(n) z^(-n)                      [5.a]
       dz          dz n=-inf

                   inf          d
                =  sum   x(n) --- z^{-n}                      [5.b]
                  n=-inf       dz

                     inf
                =  - sum   x(n)n*z^(-n-1)                     [5.c]
                    n=-inf

                
                = -z^{-1} Z{nx(n)}                            [5.d]

                = -Z{nx(n-1)}                                 [5.e]

  
------------------------------------------------------------------------

     dDFT{x(n)}       d   inf
    ------------  = ---   sum   x(n) exp(-jwn)                [6.a]
        dw           dw n=-inf

                      inf          d
                  =   sum   x(n) --- exp(-jwn)                [6.b]    
                    n=-inf        dw

                      inf
                  =   sum  -jwx(n) exp(-jwn)                  [6.c]
                    n=-inf

                  =  -jw DFT{x(n)}                            [6.d]

========================================================================

If you compare the derivative of the ZT (Eq. [5.d]) and the DFT 
(Eq. [6.d]) you find that a "z" faxtor in the ZT that has no equivalent 
"exp" factor in the DFT. I interpret that, by [3.c], as an additional 
time shift of the ZT that is not needed in the DFT.

Now, in practical cases we work with finite-length sequences, such that 
the DFT must be regardes as a zonstrained ZT. Which means that the 
practical differentiation [5.d]-[5.e], in my view, must take presedence 
over the formal [6.d].

I can't find any flaws in this argument, and I can't get a diff scheme 
based on [6.d] to work. Now, that does not mean I am right, but it 
indicates that there's more to this than apparent at first glance.

Rune

In continuous time we c an get to the Fourier Transform from the Laplace
Transform by substituting s=jw.
Similalry in discrete time we should be able to do likewise from the
(bilateral) z-TF to the DTFT by substituting z=exp(jwT). That's my
understanding.

Tom


Fred Marshall wrote:

> "Rune Allnor" <a...@tele.ntnu.no> wrote in message
> news:f...@posting.google.com...
> > M...@21cn.com (Hakuna M. C.) wrote in message
> news:<3...@posting.google.com>...
> > > Hi all,
> > >   I am using a frequency operator to act as a differentiator like
> > >
> > >   i*w <==> d/dt
> > > here w is the frequency defined in fourier space.
> >
> > Almost. What you state is valid for continuous functions. With matlab
> > (and any other FFT routine) you work with discrete sequences. You need
> > to use the differentiation theorem of the z transform. Very crudely,
> > you get something like
> >
> >    "dx(n-1)/dt" = 2*real(IFFT(k*X(k))), k=0,...,N/2
> >
> > where dx(n)/dt is a "quasi differential" of the continuous x(t).
> > Note the circular time shift!
> >
> > DISCLAIMER: I have just found out about these things, as I'm working
> > with a problem that involves carrying "quasi differentials" of discrete
> > numbers through the DFT/IDFT. There may be (and probably is) a few
> > wrong details here, but the difference between the Fourier transform
> > and z transform is cruical.
>
> Rune,
>
> What, would you say is the crucial difference?  My notion (without thinking
> about it) is that the Discrete Fourier Transform computed to yield a
> discrete spectral transform is evaluated on a unit circle and the Z
> Transform need not be but usually is anyway??
>
> Or is it that the z transform is analogous to the Laplace transform and
> systems need to be causal??
>
> Fred