proper fft of interpolated signal

[yefj]

Hello, i have a signal 3.5*10^-5 sec long as shown bellow.
It was unevenly sampled so i did cubic spine interpolation.
my signal as shown in time domain about 0.9V and 1Mhz frequency.

I want to create an accurate fft picture with Fs=10Mhz So Ts=10^-7 sec
N=signal 3.5*10^-5/10^-7=350
I want an accurate FFT with 1Hz bin 350 samples is too little.(spectral leakage)
what could be done to perform a good fft in this case?
Thanks.

from scipy.fftpack import fft
#import plotly
#import chart_studio.plotly as py
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.widgets import Cursor
from gekko import GEKKO
#%matplotlib qt
new_x=np.arange(0,3.5e-05,1/(1e7))
dataset_fft=pd.read_table("sinus_1mhz.txt")
array_fft=dataset_fft.values
m=GEKKO()
m.x=m.Param(new_x)
m.y=m.Var()
m.cspline(m.x,m.y,array_fft[:,0],array_fft[:,1])
m.options.IMODE=2
m.solve(disp=False)


fig=plt.figure()


ax=fig.subplots()
ax.grid()
cursor=Cursor(ax, horizOn=True,vertOn=True,useblit=True,color='r',linewidth =1)
#plt.plot(array_fft[:,0],array_fft[:,1])
 
plt.plot(m.x,m.y)
plt.xlabel("time")
plt.ylabel("amp")
plt.show()

 

[KlausST]

Hi,

so your sampling window for the signal is 35us ... and you want 1Hz resolution.

it´s that same as you track the route of a car for 1s and with this information you want to "calculate" it´s route for 8 hours.
Good luck.

Klaus

 

[yefj]

Hello Klauss, i understand my error.
I will take a interpolated signal of 1 second and resample it with F=10^7Hz
then i will have 1Hz for every bin.
but there is a law that the number of fft bins needs to be power of two.
the closest power of two to my 10Mhz top frequency is 2^24 which is 16 million bins.
do i have to follow this "power of two" law in my FFT?
Thanks.

 

[FvM]

I wonder what the final target of your signal processing effort is? Apparently the input signal is monofrequent. FFT isn't necessarily the best way for a precise frequency measurement.

 

[KlausST]

Hi,

This post is still valid: https://www.edaboard.com/threads/performing-high-precision-fft-in-pyhon.406503/post-1751575
But instead of reacting to it .. you decided to start a new thread.

do i have to follow this "power of two" law in my FFT?

Do you want us to read the documentation for you? I don't even know what exact FFT library you use.

Klaus

 

[yefj]

Hello Klauss , yes you are absolutely correct.you even said that using 10 million bins is too much.
that you used at most 45K for any signal.
so i interpolated a signal making it "continues" .
what are the general guidelines (general strategy) i will follow to get an accurate spectral picture.

 

[KlausST]

Hi,

I can only repeat what I´ve written: decide your true requirements.
And - for sure - if you want us to help you, you should tell us your decisions ... and why.

As also already written: I can´t see how a "1Hz" resolution on a 10MHz sampled signal should make sense.
It means 0.1ppm resolution. What hardware and what physics can reproduce that exact?
Any noise, any temperatre drift, any calculation error (even the resolution of a "float") may come to the limits.

Don´t get me wrong: from all your given information I call it nonsense.
If you have more information that make such a resolution "useful", then tell us.

--> it makes no sense to give the speed of a car with 7 digits of resolution:
like: 12.34567 km/h while the last digit means 1cm/hour. Even a snail may be 5000 times faster!

Klaus

 

[yefj]

Hello Kluass , I have built an amplifier in LTSPICE.
My goal is to see the compression of the amplifier in the frequency domain(before and after)
As we can see in the uncompressed image, our output is 1MHZ and DC component.
LTSPICE samples are not time step even so in order to make a FFT i am doing a cubic spline to get a "continues signal".
I want to see a good fft picture from DC till 10Mhz.(so we will have a mirror at 5MHZ)
In the compressed signal we have a high order tones rising.

how long should i get my "continues" interpolated signal?
how many fft bins should i use to get a good representation of my compressed and non compressed signal?
what do you reccomend?
Thanks.

not compressed:

compressed:

 

[KlausST]

Hi,

if you have a 1MHz input, low distortion, constant frequency, constant amplitude (otherwise I guess a validation makes no sense).

on a 10MHz sampling rate all you will see is:
* modified DC (bin 0)
* modified fundamental (1MHz, bin 1.000.000)
* modified overtones: 2MHz, 3MHz, 4MHz (bins: 2.000.000, 3.000.000, 4.000.000)
That´s all. No other bins are of interest.
A distortion of an 1MHz sine will cause integer multiples, so 2.000.000 Hz, (not: 1.999.999 Hz or 2.000.001 Hz).

So these are exactly 5 bins out of 5 million!!

All the other bins contain: voltage noise, frequency noise (jitter), quantisation noise, mathematical errors, your interpolation errors, maybe regulation loop operation and other artefacts.

examples:
* if you see a significant signal at 123.456 Hz there will be something wrong. It´s not caused by the input signal at all.
* If you see 1.999.923 Hz instead of 2.000.000 Hz then there will be a different problem
* If you see 1.999.923 Hz additionally to the 2.000.000 Hz then there may be a regulation loop problem (77kHz sideband caused by modulation)

So for general distortions you just need 5 bins. You don´t need a crazy "1Hz frequency resolution". Theretically you could go with a 1MHz frequency resolution.

In real world you will see "signal" in all the other 4.999.995 bins. But I guess you can gain about no information from this.
A lot of effort, but mainly useless.

Even worse:
Because of the "interpolation" you create a lot of noise (signals in other bins, as well as phase and amplitude erros in your wanted bins). Are you able to differentiate where the signals come from?
And even more worse:
If you don´t use a sampling frequency that is an exact integer number of your test signal, then you get chaotic results.
In this case one usually tries to "improve" a little by applying a windowing function (mathematical) on the FFT´s input data. But then every single frequency looks like "spread spectrum- ed" --> it will not use a single frequency bin, but a lot of neighbouring bins.

You need a lot of experience to interprete the results of the FFT. Indeed you need a lot of experience to get useful results from the FFT in frst hand.
It´s like trying a dog to control a F1 car. A F1 car can be a lot of fun, but it first needs to understand some theory and a lot of experience. Experience gained step by step.

**
I don´t want to discourage you, in no way. An FFT is an ingenious instrumet to analyze things. But you need to understand
* how to prepare the input data properly
* what to expect from the results
* how to get the results
* how to interprete the results.

Let´s imagine I want to go your way. I guess it takes me weeks of of frustration .. because the results will not be as ideal as I expect them to be. In the end I will give up. ... or lower my expectations and live with more realistic data.
More realistic, that could be obtained by much less effort within hours.

So if you see a benefit: go your way. I still try to support you with informations regarding FFT data and interpretation. But I´m not experienced with your FFT library.

Klaus

 

[kaz1]

Hello, i have a signal 3.5*10^-5 sec long as shown bellow.
It was unevenly sampled so i did cubic spine interpolation.
my signal as shown in time domain about 0.9V and 1Mhz frequency.

I want to create an accurate fft picture with Fs=10Mhz So Ts=10^-7 sec
N=signal 3.5*10^-5/10^-7=350
I want an accurate FFT with 1Hz bin 350 samples is too little.(spectral leakage)
what could be done to perform a good fft in this case?
Thanks.

from scipy.fftpack import fft
#import plotly
#import chart_studio.plotly as py
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.widgets import Cursor
from gekko import GEKKO
#%matplotlib qt
new_x=np.arange(0,3.5e-05,1/(1e7))
dataset_fft=pd.read_table("sinus_1mhz.txt")
array_fft=dataset_fft.values
m=GEKKO()
m.x=m.Param(new_x)
m.y=m.Var()
m.cspline(m.x,m.y,array_fft[:,0],array_fft[:,1])
m.options.IMODE=2
m.solve(disp=False)


fig=plt.figure()


ax=fig.subplots()
ax.grid()
cursor=Cursor(ax, horizOn=True,vertOn=True,useblit=True,color='r',linewidth =1)
#plt.plot(array_fft[:,0],array_fft[:,1])
 
plt.plot(m.x,m.y)
plt.xlabel("time")
plt.ylabel("amp")
plt.show()

View attachment 182123

I agree with various replies. Just adding few notes.
You can do any size FFT, need not be 2^n, depends on your tools.
You can do Goertzel approach for few bins instead of full FFT bins
You can repeat your signal, nicely without breaking it if you can't get 10 million samples, instead of post-ADC upsampling. Overall you will be cheating the sampling process of ADC, getting 350 but needing 10 million samples.

 

[yefj]

Hello , So given Top frequency of 10Mhz and we want to catch every multiple of a 1Mhz so we need 10 samples?
how can i get an acurate result on that 1Mhz bin sample?
As i see it, we say to the fft that we have only 10 bins and we will get a serios spectral leakage.
Or even when i do 1e4 samples.
i still will get spectral leakage.
how do i compensate the spectrl leakage and get an accurate result for 0 1Mhz 2MHz etc..?
Thanks.

 

[D.A.(Tony)Stewart]

If you want to measure phase noise with 1Hz bins from 1MHz you need at least 2 million cycles and similar samples.

If you want to measure harmonic distortion, use a spectrum analyzer of a DSO with <=10 ns sample time and process about 4k samples

I prefer to measure the sinusoidal error by measuring the difference peaks as a %ratio of the average when AC coupled or DC=0 This correlates well with THD.


Below I compare a common emitter with a gain of 18 to one that has negative feedback Rcb and distortion is reduced from 30% to 3% approx. and gain reduced to 10.

Vbe modulation has 2nd order effect on THD from reduced gain near cutoff and towards clipping causing odd harmonic distortion such as in your plot.

Since you had an XY question about FFT, what are ALL your design specs and challenges?

FWIW

 

[KlausST]

Hi,

What do you mean with spectral leakage?

Sampling rate of 10MSmpl/s and 10 samples:
--> you get 0Hz, 1MHz, 2MHz, 3MHz, 4MHz. That's it. 5 output bins.

Maybe you want ask what happens if there is a frequency of 4.25Hz?
--> most important answer: it violates the FFT theory ... as an FFT input is meant as a countinous signal repeating itself (after the 10 samples).
But if you "repeat" this 10 samples you will see an "unnatural" jump from last sample to first sample.
--> so you will see the highest peak at the 4MHz bin ... but because of the jump you will also get "invalid" signal in all other bins.

So if you really need a 0.25MHz resolution on a 10MSmpl/s rate you need 40 samples. But 40 is not a power of 2 number.
So either you need to use an FFT algorithm that can work with _non_ power of 2 count, or you need to use a windowing function.
But a windowing function will modify the result ... regarding sharpness of the frequency resolution.

******
For THD measurement applications you should use
* one common low jitter high frequency clock
* a DAC for "synchronous" frequency generation, (sample rate drived from above common clock)
* an ADC for "synchronous" sampling, (sample rate drived from above common clock)
Then you can guarantee that the measurement signal is synchronous to sampling (and even the power of 2 count)
You don't need a windowing function
While you get perfect, clean, sharp FFT results.
The errors (DAC aperture jitter, DAC noise ... ADC errors, .... clock jitter ....) determine the overall performance.

For me the measurement precision (=repeatability) tells me how reliable the measurement results are.
Thus I recommend you the same what I do: Check the reliability of your measurement results.

Do a real measurement: perform all that is necessary to get the result. Archive the result.
Then do a complete new masurement ... and compare the result with the other.
Reapeat this seveal times.
(Assuming your 10Msamples/s for 1 second)
Now you will see that one has a peak at 4,000,009 Hz and the other is at 3,999,994Hz.
--> so how "useful" is a 1Hz resultion when it differs -6 ... +9Hz? You simply can't rely on the "1Hz" information. (Drift thermal/voltage/ageing drift of system clock frequency or measurement frequency)

Also you will see that there it shows -7.8dB, while the other shows -8.2dB ... again an uncertainty of 0.4dB.

Also you will see that the next bins beside the "peak" bin also show a lot of signal. Caused by sample clock jitter, measurement signal jitter, mathematical interpolation problems...)

... all while you know that the overtone of a 1,000,000 signal can only be a clean integer multiple = 4,000,000 MHz.

****
Imagine you buy a dozen screws. But a weigh scale tells this are 12.0436 screws. Buy another 12 scriws and the weigh scale shows 11.996 dcrews.
So how important (informative) can the 4 decimals be? For your application. Most probably you need them just to be 12.
It depends on application, because the quality control of a screw manufacturer may gain some information from decimals ... but maybe only of 2 decimals.

Klaus

 

[yefj]

Hello Klauss,In my interpolated signal i have resampled it every 1/1e7 and i got 793 samples.
I have done zero padding to the left 100000-793 zeros, my fft plot havent changed at all after padding.
my original file shown in the link bellow.
the edges are not equation so the sudden discontinuty creates noise to all the spectrum.
My python code is shown bellow.

my main tones is DC and
I am trying to make Keyser window to improve the situation .
How to you reccomend to build the keiser taper?
Thanks.

https://numpy.org/doc/stable/reference/generated/numpy.kaiser.html

original data:

sinus_1mhz

Visit this link to preview document

files.fm

# -*- coding: utf-8 -*-
"""
Created on Thu Mar 30 13:04:11 2023

@author: Asus
"""
import math
from scipy.fftpack import fft
#import plotly
#import chart_studio.plotly as py
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.widgets import Cursor
from gekko import GEKKO
#%matplotlib qt
new_x=np.arange(0,7.89752e-05,1/(1e7))
dataset_fft=pd.read_table("sinus_1mhz.txt")

array_fft=dataset_fft.values
m=GEKKO()
m.x=m.Param(new_x)
m.y=m.Var()
m.cspline(m.x,m.y,array_fft[:,0],array_fft[:,1])
m.options.IMODE=2
m.solve(disp=False)


freq_vec=1e7*np.arange(0,1,1/(1e7))

m2=np.pad(m.y, (0,100000-len(m.y)), 'constant')

fft_y=fft(m2,len(freq_vec))



fig=plt.figure()

ax=fig.subplots()
ax.grid()
cursor=Cursor(ax, horizOn=True,vertOn=True,useblit=True,color='r',linewidth =1)
plt.plot(freq_vec,20*np.log10(abs(fft_y)/max(abs(fft_y))))
plt.xlabel("Freq")
plt.ylabel("amp_normalised")
plt.show()

 

[kaz1]

Your file shows 1993 samples (not 793).

I removed dc by subtracting mean then did fft of 2k (i.e. zeros = 2048-1993):

 

[yefj]

Hello Kaz1,the amplitude is not 90 of the sine.I think also your frequency tone is not at 1Mhz.

 

[kaz1]

I removed dc offset then did fft without descaling. FFT tells me your tone is @ ~390KHz assuming Fs is 10Msps.

 

[yefj]

Hello Kaz1, I have shown already in the results post before that my tone in fft is 1MHZ and the ampitude is far higher then the original.
I showd my fft i got a similar photo and more accurate with the code i posted in python.
But the purpose is to know what kind of windowing to use in this case?
Thanks.

 

[kaz1]

if your Fs is 10Msps then the vector I checked is wrong. I did not care about scaling. The vector I got is + 1.17~ +0.97 as below, tone @ 395KHz:

 

[D.A.(Tony)Stewart]

I removed dc offset then did fft without descaling. FFT tells me your tone is @ ~390KHz assuming Fs is 10Msps.

The raw data tells me it is 1/8 MHz and 1 MHz doing my eyeball inversion of 20 us of raw data.

Why? assuming the X samples are not evenly spaced and represent relative time.