I'm trying to get fft amplitudes to match the input time domain signals amplitude scale.
I've tried so many 'answers' to this problem - none of which seem to work for me. So I've written a program to experiment.
This takes a section of audio (from a file or - you can use additative synthesis to create a wave).
I then give the option of padding the signal with zeros (choose number and number of steps to analyse).
I don't have a site to post images of the plots I've created - but if you cut and post the program at the bottom (python 3.x) you'll quickly see them for yourselves!
A small sub-section of the program is below.
The pertinent section of code (plot_pad_stats() below) is:
for p in np.arange(start,end,step):
s1 = pad(s,int(p))
ftp = np.fft.fft(s1,norm=norma)
ft = norm_fft(s1,norm)
if real:
ftp = ftp.real
ft = ft.real
else:
ftp = np.sqrt(ftp.real**2 + ftp.imag**2)
ft = np.sqrt(ft.real**2 + ft.imag**2)
if absolute:
ftp = abs(ftp)
ft = abs(ft)
if active:
ftp = [(f+abs(f))/2 for f in ftp]
ft = [(f+abs(f))/2 for f in ft]
minft[p1] = min(ftp)
maxft[p1] = max(ftp)
meanft[p1] = np.mean(ftp)
lenft[p1] = int(p)
minnft[p1] = min(ft)
maxnft[p1] = max(ft)
meannft[p1] = np.mean(ft)
lennft[p1] = int(p)
I suspect something to do with phase - but have no idea how to get an fft amplitude scale to match the time domain signal being analysed...
Not only is the change as padding is added interesting, but the change as the signal goes from a single sine wave to more complex waves!
Any help gratefully appreciated to:- tony.fdlr@gmail.com
arrays = plot_pad_stats: just call with no args to plot 440hz signal full amplitude
use
additativeSignal(hza=[1],aa=None,ap=[0],n=88200,sr=44100,verbose=True):
'''
hza = herz array
aa = amplitude array
ap = phase array
n = number samples to return
sr = sampling frequency
verbose = sanity check
'''
to create additative synthesised signals eg
signal = additativeSignal([440,880,1320],[0.6,0.3,0.2],[0,0,0])
signal = additativeSignal([440,550,600],[0.6,0.5,0.3],[0,50,100])
then
plot_pad_stats(signal)
use
plot_signal(s,start,end) to plot the signal (or a section thereof)
Note = phase values are in samples not degrees or radians!
OR for more fun:
sawtooth wave:
f = 100
a = 0.8
ph = 0
noFqs = 20
signal = additativeSignal(
[f * i for i in range(1,noFqs+1)],
[am/i for i in range(1,noFqs+1)],
[0 for i in range(1,noFqs+1)])
plot_signal(signal[0:1000])
arrs = plot_pad_stats(signal,maxpad=1000,steps=1000,absolute=False)
Play with the other args to your heart's content!
Be aware it can take a little time to plot!
Here's the code for you to play with (python 3.x)
import numpy as np
import matplotlib.pyplot as plt
def pad(s,n):
'''
pad signal to make length n
'''
ls = len(s)
#print('Pad request: sig len: {}, new len: {}'.format(ls,n))
if n <= ls:
return s
s2 = np.zeros(n)
try:
s2[:ls] = s
except:
#print('Pad exception: sig len: {}, new len: {}'.format(ls,n))
return s
return s2
def norm_fft(s,norm=True):
'''
See plot_pad_stats below for explanation
'''
if norm:
norma = 'ortho'
else:
norma = None
n = len(s)
ft = np.fft.fft(s,norm=norma) / n
ft = ft[range(int(n / 2))]
return ft
def additativeSignal(hza=[1],aa=None,ap=[0],n=88200,sr=44100,verbose=True):
'''
hza = herz array
aa = amplitude array
ap = phase array
n = number samples to return
sr = sampling frequency
verbose = sanity check
'''
xa = np.arange(0,n)
s = np.zeros(n)
a = 0
apl = len(ap)
if apl < len(hza):
apl = np.zeros(len(hza))
if not (isinstance(aa,(np.ndarray,list))):
aa = np.ones(len(hza))
if not (isinstance(ap,(np.ndarray,list))):
aa = np.ones(len(hza))
for hz in hza:
amp=aa[a]
phase = ap[a]
if verbose:
print("freq %d, amp = %f, phase: %d" % (hz,amp,phase))
if hz != 0:
s += (np.sin(2*np.pi*(xa+phase)/(sr/hz))*amp)
a += 1
return s
def plot_pad_stats(s=None,
sr=44100,
minpad = 0,
maxpad = 300,
steps = 300,
real = True,
norm = None,
absolute = True,
active = False
):
'''
s = signal (audio etc)
sr = sample rate
minpad - minimum number of zeroed samples to append
maxpad - maximum number of zeroed samples to append
steps - number of ffts to calculate
real - use only real values otherwise take sqrt(real**2 + imag**2)
norm - np.fft.fft norm arg
absolute - make fft results absolute (ie >0)
active - standard activation function - any negative values changed to zero
'''
if norm:
norma = 'ortho'
else:
norma = None
minpad = int(minpad)
# Create 1 sec 4450 herz full volume audio signal at 44100 hz sampling frequency
if type(s) == type(None):
s = additativeSignal([440],[1],[0],44100,44100,True)
start = len(s) + minpad
end = start + maxpad
prange = maxpad
#Arrays of fft results:
# length fft
lenft = [0 for p in range(steps)]
# max fft val
maxft = [0 for p in range(steps)]
# min fft val
minft = [0 for p in range(steps)]
# mean fft val
meanft = [0 for p in range(steps)]
# same arrays using the above normalised function
lennft = [0 for p in range(steps)]
maxnft = [0 for p in range(steps)]
minnft = [0 for p in range(steps)]
meannft = [0 for p in range(steps)]
p1 = 0
step = prange/steps
print('fft len sig {} step {}'.format(len(s),step))
for p in np.arange(start,end,step):
#print('fft {} len sig {} padded len {}'.format(p1,len(s),int(p)))
s1 = pad(s,int(p))
ftp = np.fft.fft(s1,norm=norma)
ft = norm_fft(s1,norm)
if real:
ftp = ftp.real
ft = ft.real
else:
ftp = np.sqrt(ftp.real**2 + ftp.imag**2)
ft = np.sqrt(ft.real**2 + ft.imag**2)
if absolute:
ftp = abs(ftp)
ft = abs(ft)
if active:
ftp = [(f+abs(f))/2 for f in ftp]
ft = [(f+abs(f))/2 for f in ft]
minft[p1] = min(ftp)
maxft[p1] = max(ftp)
meanft[p1] = np.mean(ftp)
lenft[p1] = int(p)
minnft[p1] = min(ft)
maxnft[p1] = max(ft)
meannft[p1] = np.mean(ft)
lennft[p1] = int(p)
p1 += 1
#:22.8f
plt.subplot(211)
plt.title("FFT stats. Sig Len {} with {} - {} appended zeros".format(len(s),minpad,maxpad+minpad))
plt.xlabel("signal + pad length {} - {}".format(start,end))
plt.ylabel("max {:f} - {:f} and\nmin {:f} - {:f}\nand mean {:f} - {:f}\nffta values (abs/real)".format(min(maxft),max(maxft),min(minft),max(minft),min(meanft),max(meanft)))
plt.plot(lenft,maxft,'bo')
plt.plot(lenft,minft,'r+')
plt.plot(lenft,meanft,'g*')
plt.subplot(212)
plt.title("\nNormalised FFT stats. Sig Len {} with {} - {} appended zeros".format(len(s),minpad,maxpad+minpad))
plt.xlabel("signal + pad length\n{} - {}".format(start,end))
plt.ylabel("max {:22.8f} - {:22.8f} and\nmin {:22.8f} - {:22.8f}\nand mean {:22.8f} - {:22.8f}\nffta values (abs/real)".format(min(maxnft),max(maxnft),min(minnft),max(minnft),min(meannft),max(meannft)))
plt.plot(lennft,maxnft,'bo')
plt.plot(lennft,minnft,'r+')
plt.plot(lennft,meannft,'g*')
plt.show()
return maxft,minft,lenft
def plot_signal(s,start=0,end=None):
plt.title("Signal")
plt.xlabel("samples")
plt.ylabel("amplitude")
if end == None:
plt.plot(s)
else:
plt.plot(s[start:end])
plt.show()
Regards
Tony
