Получить амплитуды БПФ для соответствия шкале амплитуд входных сигналов во временной области (numpy.fft) - PullRequest
0 голосов
/ 12 ноября 2019

Я пытаюсь получить амплитуды БПФ в соответствии со шкалой амплитуды входных сигналов во временной области.

Я пробовал так много «ответов» на эту проблему - ни один из которых мне не подходит. Итак, я написал программу для экспериментов.

Для этого требуется фрагмент аудио (из файла или - вы можете использовать дополнительный синтез для создания волны).

Затем я даю вариантзаполнения сигнала нулями (выберите количество и количество шагов для анализа).

У меня нет сайта для публикации изображений созданных мною графиков, но если вы вырезаете и публикуете программу ввнизу (python 3.x) вы сами быстро их увидите!

Небольшой подраздел программы находится ниже.

Соответствующий раздел кода (plot_pad_stats () ниже):

    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)

Я подозреваю, что что-то связано с фазой, но не знаю, как получить шкалу амплитуды БПФ, соответствующую анализируемому сигналу во временной области ...

Не толькоЭто изменение, так как дополнение добавляется интересно, но изменение по мере того, как сигнал переходит от одной синусоиды к более сложным волнам!

Любая помощь с благодарностью: - tony.fdlr@gmail.com

arrays = plot_pad_stats: просто позвонитенет аргументов для построения полной амплитуды сигнала 440 Гц

используйте

   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])

затем

    plot_pad_stats(signal)

используйте

    plot_signal(s,start,end) to plot the signal (or a section thereof)

Примечание = значения фазы приведены в образцах, а не в градусах или радианах!

ИЛИ для большего удовольствия:

пилообразная волна:

        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)

Играйте с другими аргументами, которые вам по душеcontent!

Помните, что заговор может занять некоторое время!

Вот код, с которым вы можете играть (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()

С уважением

Тони

1 Ответ

0 голосов
/ 12 ноября 2019

Если я вас правильно понял, вы хотите увидеть правильную амплитуду временного сигнала в спектральной области. В приведенном ниже примере я добавляю две синусоидальные волны амплитуд 1 и 0,5. В спектре справа вы можете видеть, что амплитуды выглядят вдвое, что связано со сложным БПФ, разделяющим мощность между отрицательной и положительной частотами.

import numpy as np
import matplotlib.pyplot as p
%matplotlib inline

T=3 # secs
d=0.04 # secs
n=int(T/d)
print(n)
t=np.arange(0,T,d)  
fr=1 # Hz
y1= np.sin(2*np.pi*fr*t) 
y2= 1/2*np.sin(2*np.pi*3*fr*t+0.5) 
y=y1+y2 
f=np.fft.fftshift(np.fft.fft(y))
freq=np.fft.fftshift(np.fft.fftfreq(n,d))
p.figure(figsize=(12,5))
p.subplot(121)
p.plot(t,y1,'.-',color='red', lw=0.5, ms=1)  
p.plot(t,y2,'.-',color='blue', lw=0.5,ms=1) 
p.plot(t,y,'.-',color='green', lw=4, ms=4, alpha=0.3) 
p.xlabel('time (sec)')
p.ylabel('amplitude (Volt)')
p.subplot(122)
p.plot(freq,np.abs(f)/n)
p.xlabel(' freq (Hz)')
p.ylabel('amplitude (Volt)');

enter image description here

...