Author Topic: FFT and Signals  (Read 2605 times)

0 Members and 1 Guest are viewing this topic.

Offline LangerTopic starter

  • Newbie
  • Posts: 6
  • Country: de
FFT and Signals
« on: March 15, 2019, 07:16:21 pm »
Hello,

I have a question regarding DFT/Fast Fourier Transformation and signal processing.

I want to write a program (matlab in this case) which can calculate the reaction (current) on an arbitrary waveform (voltage).

My idea was to take the arbitrary waveform, do FFT, devide the FFT coefficients with the impedance at the desired frequency: I(f)=U(f)/Z(f), and then do an inverse FFT I(f)->I(t).

This is working fine for harmonic waveforms. But when I try to do this with a square wave it fails at the edges.

The attached picture shows U(t), I(t) calculated with matlab and the same thing with ltspice. The rising current at the edge seems to be wrong comparing the result with ltspice.


Thank you in advance for every hint and solution.

matlab code:

======================
clear all;
clc;

%define measurement
Fs = 10000;            % Sampling frequency                   
T = 1/Fs;             % Sampling period       
L =Fs;             % Length of signal
t = T*(0:L-1);        % Time vector
f = Fs*(1:L-1)/L;

%define input waveform U(t)
U = 4*(1+square(10*t))/2;
%U0=2;
%U = U0*exp(i*2*pi*300*t)+U0*exp(i*2*pi*100*t)


%do the fft of U(t)->U(f)
Ufft = fft(U);

%plot the FFT of U(t)
fig_4=figure(4);
ax_4=axes;
semilogy(abs(Ufft));

%plot U(t)
fig_1=figure(1);
ax_1=axes;
plot(ax_1,t,real(U))
title('Input U(t)')
fig_1.Position = [0,343,560,420];

%define a frequency vector with rising and falling frequencies looking like
%this: /\ because of the mirror coefficients.
freq=Fs*(0:L/2-1)/L;
freq=horzcat(freq,fliplr(freq));

%define the impedance (here only C and R in series)
L=0;
Cap=100e-6;
R=1;

res=1/(2*pi*sqrt(L*Cap))    %calculate resonance frequency
Xl=i*2*pi*L.*freq;
Xc=-i./(2*pi*Cap*freq);
Z = R+Xl+Xc;                %calculate impedance of the system


fig_2=figure(2);
ax_2=axes;
semilogy(freq,abs(Z))
title('Impedance Spectrum of Z(f)')
fig_2.Position = [560,343,560,420];


Ifft=Ufft./Z;       %devide the Fourier coefficients by the impedance Z(f)
fig_3=figure(3);
ax_3=axes;
current=ifft(Ifft);
plot(ax_3,t,real(current))
fig_3.Position = [1120,343,560,420];
title('Output I(t)')
max(current)

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

Offline T3sl4co1l

  • Super Contributor
  • ***
  • Posts: 22436
  • Country: us
  • Expert, Analog Electronics, PCB Layout, EMC
    • Seven Transistor Labs
Re: FFT and Signals
« Reply #1 on: March 15, 2019, 07:55:07 pm »
Without reviewing the code -- looks like it's just missing a U(t) (unit step function).  Fourier is for all time so you need to check your conditions as far as events (zeroing negative time?) and such go.

Tim
Seven Transistor Labs, LLC
Electronic design, from concept to prototype.
Bringing a project to life?  Send me a message!
 

Online radiolistener

  • Super Contributor
  • ***
  • Posts: 4135
  • Country: 00
Re: FFT and Signals
« Reply #2 on: March 15, 2019, 08:15:03 pm »
At a glance it's because you're scaling frequency magnitudes with variable coefficient Z.
May be something wrong with Z?

PS: interesting approach for current simulation  :) But I'm not sure, is it correct to do it in such way?
« Last Edit: March 15, 2019, 08:39:55 pm by radiolistener »
 

Offline LangerTopic starter

  • Newbie
  • Posts: 6
  • Country: de
Re: FFT and Signals
« Reply #3 on: March 15, 2019, 10:19:10 pm »
Without reviewing the code -- looks like it's just missing a U(t) (unit step function).  Fourier is for all time so you need to check your conditions as far as events (zeroing negative time?) and such go.

Can you specify this? When i look at FFT I see the amplitudes twice ( matlab shows them only in positiv frequency direction and not negative). So i should delete the doubled/negative frequencies?
 

Offline LangerTopic starter

  • Newbie
  • Posts: 6
  • Country: de
Re: FFT and Signals
« Reply #4 on: March 15, 2019, 10:24:42 pm »
At a glance it's because you're scaling frequency magnitudes with variable coefficient Z.
May be something wrong with Z?

PS: interesting approach for current simulation  :) But I'm not sure, is it correct to do it in such way?

I'm pretty sure that Z is fine.
Deivding the freq amplitude U by Z sounds reasonable for me. But if I'm wrong please correct me....
And it works for harmonic signals, so it seems correct.
 

Offline unitedatoms

  • Frequent Contributor
  • **
  • !
  • Posts: 324
  • Country: us
Re: FFT and Signals
« Reply #5 on: March 15, 2019, 10:35:09 pm »
The Fourier does not care about causality. It also does not care about finite span of observed signal. It thinks that signal was repeated over eternity. So to perform closer to intuitive physics, try to do windowing, and also try to shift phase of pulse to some odd angle. And ignore reconstructed "predicted past", ignore negative frequencies etc.
Interested in all design related projects no matter how simple, or complicated, slow going or fast, failures or successes
 

Offline dmills

  • Super Contributor
  • ***
  • Posts: 2093
  • Country: gb
Re: FFT and Signals
« Reply #6 on: March 15, 2019, 10:59:39 pm »
Looking at the code, it looks like your 'square' wave is at Fs/10, which will alias like a bastard (It only fits ONE harmonic, and the 5th aliases to DC)..

Regards, Dan.
 

Offline bson

  • Supporter
  • ****
  • Posts: 2497
  • Country: us
Re: FFT and Signals
« Reply #7 on: March 15, 2019, 11:30:22 pm »
Deivding the freq amplitude U by Z sounds reasonable for me
Did I miss something important here?  Multiplication in the time domain is equivalent to convolution in the frequency domain and vice versa. In the frequency domain,  I=U*1/Z becomes F(I) = F(U*1/Z) = F(U)**F(1/Z) where ** is the convolution operator.   So instead of dividing Ufft with the impedance, shouldn't you convolve with its reciprocal?  (Of course, is Z is constant rather than Z(w), then F(1/Z) is constant, but then you also have no reactive current.)
« Last Edit: March 15, 2019, 11:32:04 pm by bson »
 

Offline T3sl4co1l

  • Super Contributor
  • ***
  • Posts: 22436
  • Country: us
  • Expert, Analog Electronics, PCB Layout, EMC
    • Seven Transistor Labs
Re: FFT and Signals
« Reply #8 on: March 16, 2019, 02:08:17 am »
Without reviewing the code -- looks like it's just missing a U(t) (unit step function).  Fourier is for all time so you need to check your conditions as far as events (zeroing negative time?) and such go.

Can you specify this? When i look at FFT I see the amplitudes twice ( matlab shows them only in positiv frequency direction and not negative). So i should delete the doubled/negative frequencies?

Specifically, frequency domain needs to be complex conjugate (which will look equal amplitude if you're only looking at absolute value) to get an all-real output; but I don't remember what you need to get a unique negative-zero-valued output.  In any case, if the math is correct, then you can ignore t < 0.

Tim
Seven Transistor Labs, LLC
Electronic design, from concept to prototype.
Bringing a project to life?  Send me a message!
 

Offline rhb

  • Super Contributor
  • ***
  • Posts: 3516
  • Country: us
Re: FFT and Signals
« Reply #9 on: March 16, 2019, 02:23:16 am »
1)  The discrete (sampled) Fourier transform is periodic over [-pi:pi) or [0:2*pi) depending upon your conventions.  The latter is what you actually get out of an FFT.   The negative frequncies appear on the right.

2) To make a signal causal, i.e. zero before T0, the imaginary part needs to be the Hilbert transform of the real part.  In the frequency domain the real part is even and the imaginary is odd.  The usual term is the transform is hermitian.  The transform of a hermitian function is pure real and causal.

Compute the FFT of a zero phase (symmetric about T0) function.  The imaginary part will be uniformly zero.  Multiply the negative frequencies by -1, the positive frequencies by 1 and set zero frequency to zero and store that in the imaginary part of the transform.  then inverse transform.  The result will be causal and minimum phase.
 

Offline T3sl4co1l

  • Super Contributor
  • ***
  • Posts: 22436
  • Country: us
  • Expert, Analog Electronics, PCB Layout, EMC
    • Seven Transistor Labs
Re: FFT and Signals
« Reply #10 on: March 16, 2019, 02:43:16 am »
That's right, Hermitian is what I was thinking of.

Tim
Seven Transistor Labs, LLC
Electronic design, from concept to prototype.
Bringing a project to life?  Send me a message!
 

Offline LangerTopic starter

  • Newbie
  • Posts: 6
  • Country: de
Re: FFT and Signals
« Reply #11 on: March 16, 2019, 01:18:02 pm »
1)  The discrete (sampled) Fourier transform is periodic over [-pi:pi) or [0:2*pi) depending upon your conventions.  The latter is what you actually get out of an FFT.   The negative frequncies appear on the right.

2) To make a signal causal, i.e. zero before T0, the imaginary part needs to be the Hilbert transform of the real part.  In the frequency domain the real part is even and the imaginary is odd.  The usual term is the transform is hermitian.  The transform of a hermitian function is pure real and causal.

Compute the FFT of a zero phase (symmetric about T0) function.  The imaginary part will be uniformly zero.  Multiply the negative frequencies by -1, the positive frequencies by 1 and set zero frequency to zero and store that in the imaginary part of the transform.  then inverse transform.  The result will be causal and minimum phase.

2) Let me summarize what I understood:
- Make the input Signal U(t) symmetric to t=0; (the imag part is zero with some tricks in matlab (I double checked that));
- multiply pos. by 1 and neg. freq by -1 (okay.....I could do it)
- set zero freq to zero is unnecessary because it already should be zero (but I could do it)

But I don't understand your next point: 'store that in the imaginary part of the transform' store what in the imag part?!

And the other question: Is it wright to devide the FFT signal by Z and then transform it back (this might be the most important question)?!

 

Offline rhb

  • Super Contributor
  • ***
  • Posts: 3516
  • Country: us
Re: FFT and Signals
« Reply #12 on: March 16, 2019, 02:03:26 pm »
Given f(t) for a real function which is symmetric about T0 (e.g. the autocorrelation) compute the transform, F(w).  That will  be a pure real complex series.  Compute the Hilbert transform of the real part as described and store that in the imaginary part of F(w).  Then back transform.  The result will be pure real and zero for all times before T0.

In the general case F(0) is not zero.  So setting the imaginary part of F(0) to zero is important.

It would probably help to review the mathematics of the analytic signal.  What is being done above is the mathematical equivalent of creating the quadrature part of an I/Q series.  The Hilbert transform is a 90 degree phase shift operator.

BTW don't feel badly about your confusion on this.  It is often badly described.  When I was at Austin working under one of the members of Norbert Wiener's Geophysical analysis Group, none of us could figure out how to make a minimum phase wavelet this way.  So we just applied a Wiener inverse twice.  I only figured out how to do it via the Hilbert transform in the last year.  After 30 years of industry DSP experience!   It rarely came up and if it did I was in too much of a hurry to finish the work, so I just used the Wiener inverse dodge.

I can't answer about the impedance.  I'm a geophysicist by training, so I'm just getting started doing this stuff in the EM field.  I was always working with elastic wave propagation and we were only concerned with the reflection coefficients.  There are 3 of those one compressional and two shear, all of which vary with angle of incidence.  That's for isotropic, plane layered media.  Reality is anisotropic, so the coefficients also vary with azimuth in the general case which is what one usually has.  And it's not really plane layers.
 

Offline LangerTopic starter

  • Newbie
  • Posts: 6
  • Country: de
Re: FFT and Signals
« Reply #13 on: March 16, 2019, 03:01:39 pm »
I found some information about the Hilbert Transformation: https://www.gaussianwaves.com/2017/04/analytic-signal-hilbert-transform-and-fft/

In general my approach seems to work with harmonic functions, but fails when going to square wave function.

Even if I do the hilbert transformation, I'm not sure if Ufft/Z is then still applicable.

So I will try a new approach which is called convolution: https://www.dspguide.com/ch13/2.htm
This approach seems to fit my needs.

But thanks everybody for his/her hints!

I will post some code if I get it running...

 

Offline rhb

  • Super Contributor
  • ***
  • Posts: 3516
  • Country: us
Re: FFT and Signals
« Reply #14 on: March 16, 2019, 03:19:53 pm »
Square waves aren't a problem except that you have to deal with the zeros of sin(x)/x.  That's usually done by making them a small percentage of the modulus of the peak value with the correct phase.
 

Offline bson

  • Supporter
  • ****
  • Posts: 2497
  • Country: us
Re: FFT and Signals
« Reply #15 on: March 17, 2019, 12:41:27 am »
I posited:
(Of course, is Z is constant rather than Z(w), then F(1/Z) is constant, but then you also have no reactive current.)
On second thought this is not true at all.  Convolution with a fixed, DC kernel (1/Z) still requires convolution.
 

Offline LangerTopic starter

  • Newbie
  • Posts: 6
  • Country: de
Re: FFT and Signals
« Reply #16 on: April 04, 2019, 05:33:52 pm »
I finally made a working matlab-code.

This can calculate the response of an arbitrary system by using "lsim" code from matlab. Description is here: https://www.mathworks.com/help/control/ref/lsim.html
It is about convolution, which was mentioned before.

One book I found which covers this topic:
https://books.google.de/books?id=wnpTDwAAQBAJ&printsec=frontcover&hl=de

but its quite expensive...

Thanks everyone for his contribution!

matlab code:
=======

clear all;
close all;
clc;

Cap=1e-6;
R1=100;
L=0.01;


s = tf('s');
Z1=R1;
Z2=1/(s*Cap);
Z3=s*L;
%Voltage across R1:
Hs=Z1/(Z1+Z2+Z3);


t = 0:0.00001:0.01;
u = 0.5*(square(t*2*pi*500)+1);
y=lsim(Hs,u,t);
%Calculate current for R1:
y=y/R1;

fig_1=figure(1);
ax_1=axes;
plot(ax_1,t,u)
title('Voltage input')
fig_1.Position = [100,343,560,420];

hold 'on'
yyaxis right
plot(t,y)


 

Offline pwlps

  • Frequent Contributor
  • **
  • Posts: 372
  • Country: fr
Re: FFT and Signals
« Reply #17 on: April 05, 2019, 01:25:08 pm »
Deivding the freq amplitude U by Z sounds reasonable for me
Did I miss something important here?  Multiplication in the time domain is equivalent to convolution in the frequency domain and vice versa. In the frequency domain,  I=U*1/Z becomes F(I) = F(U*1/Z) = F(U)**F(1/Z) where ** is the convolution operator.   So instead of dividing Ufft with the impedance, shouldn't you convolve with its reciprocal?  (Of course, is Z is constant rather than Z(w), then F(1/Z) is constant, but then you also have no reactive current.)

What you missed is that it is the opposite: here we have convolution in the time domain which is equivalent to multiplication in the frequency domain.  I=U*1/Z  stands for I(f)=U(f)*1/Z(f). 
 

Offline pwlps

  • Frequent Contributor
  • **
  • Posts: 372
  • Country: fr
Re: FFT and Signals
« Reply #18 on: April 05, 2019, 03:21:04 pm »
1)  The discrete (sampled) Fourier transform is periodic over [-pi:pi) or [0:2*pi) depending upon your conventions.  The latter is what you actually get out of an FFT.   The negative frequncies appear on the right.

2) To make a signal causal, i.e. zero before T0, the imaginary part needs to be the Hilbert transform of the real part.  In the frequency domain the real part is even and the imaginary is odd.  The usual term is the transform is hermitian.  The transform of a hermitian function is pure real and causal.

Compute the FFT of a zero phase (symmetric about T0) function.  The imaginary part will be uniformly zero.  Multiply the negative frequencies by -1, the positive frequencies by 1 and set zero frequency to zero and store that in the imaginary part of the transform.  then inverse transform.  The result will be causal and minimum phase.


It is not the signal that should be causal but the transfer function, which is 1/Z here.
Looking at the code it is obviously hermitian :

Xl=i*2*pi*L.*freq;
Xc=-i./(2*pi*Cap*freq);
Z = R+Xl+Xc; 

therefore there is no problem related to causality here. 

One thing to check is the behavior at zero frequency (Z will diverge for freq=0), I actually don't know how Matlab handles infinite values, but even if it handles them there will be a huge numerical artifact unless you somehow cure the singularity.

I don't have Matlab but I checked with an equivalent code under Igor Pro, it works if I prevent Z from going to infinity at zero frequency, for example adding a big resistance in parallel with C. I post the results and the corresponding code in the pictures below. 


More generally,  calculations using the FFT should give the correct result if the boundary conditions are chosen to not to generate artifacts.  As unitedatoms pointed out:

The Fourier does not care about causality. It also does not care about finite span of observed signal. It thinks that signal was repeated over eternity. So to perform closer to intuitive physics, try to do windowing, and also try to shift phase of pulse to some odd angle. And ignore reconstructed "predicted past", ignore negative frequencies etc.

therefore it is important to have in mind that if the signal is not the same at start and at end then  you create an artificial discontinuity at t=0.  The same effect of periodicity appears in the simulation I show. The function U(t) is condidered as periodic and here the period of U(t) is 10, U(0)=1 and U returns to 0 at t=3*pi=9.42 which is equivalent to t=9.42-10=-0.58.  Therefore at t=0 there will be a signal reminiscent from the transient at t=-0.58, as seen in the picture of I(t).
« Last Edit: April 05, 2019, 03:39:38 pm by pwlps »
 

Offline rhb

  • Super Contributor
  • ***
  • Posts: 3516
  • Country: us
Re: FFT and Signals
« Reply #19 on: April 05, 2019, 04:38:23 pm »
 If x(0) != x(n) then things get more complicated.  In filtering operation subtracting a linear function, filtering and adding back the linear function has been my normal practice doing fiddly stuff.  For routine work a cosine taper edge is pretty common.

The most common mistake I have seen is failing to pad the series so that wraparound occurs.
 


Share me

Digg  Facebook  SlashDot  Delicious  Technorati  Twitter  Google  Yahoo
Smf