问题描述

python 这样获取信号的频率是可以的，但是不知道振幅和相位为什么不对啊？

 #!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
@author: akon
"""
from __future__ import division
import matplotlib.pyplot as pl
import numpy as np
from scipy.signal import butter, lfilter
import wave
import struct
def band_pass(freq,sig,Fs,N):
   fft_sig = np.fft.rfft(sig,N)/N
   freqs=np.fft.fftfreq(len(fft_sig))
   nyq=0.5*Fs
   j=0
   signature=[None]*len(freq)
   for i in freq:
    print i
    low=(i-10)/nyq
    high=(i+30)/nyq
    b,a=butter(6,[low,high],btype='band')
    signature[j]=lfilter(b,a,sig)
    j=j+1
   #print(signature)

   return  signature

def freq_amp(sig,N,Fs):
    fft_sig = np.fft.rfft(sig)/len(sig)
    freqs=np.fft.fftfreq(len(fft_sig))
    amp=2*np.abs(fft_sig)

    fft_amp=np.where(amp>0.3,1,0)
    print fft_amp
    #idx = [None] * N
    j=0
    m=0
    n=0
    q=0
    indx=[None]*len(fft_sig)
    for i in range(len(fft_sig)):###获取振幅为1的下标数组，并且唯一
        if((fft_amp[i]==1) and(fft_amp[i-1]-fft_amp[i]<0)):
                indx[j]=i
                j=j+1
    freq=[None]*j
    amplitude=[None]*j
    phase=[None]*j
    for k in range(j):###根据频谱振幅的下标获取频率，振幅
            freq[m]=indx[k]*Fs/N
        #freq=abs(freqs[indx[k]]*framerate)
        print indx[k]
            m=m+1
            amplitude[n]=2*np.abs(fft_sig[indx[k]])
        n=n+1
        phase[q]=np.rad2deg(np.angle(fft_sig[indx[k]]))
        q=q+1
    print freq,amplitude,phase
        return freq,amplitude,phase,indx,j

#def open(path):
#    wav_file= wave.open(path, 'rb')
#    Fs = wav_file.getframerate()
#    params = wav_file.getparams()
#    nchannels, sampwidth, framerate, nframes = params[:4]
#    data=wav_file.readframes(nframes)
#    wav_file.close()
#    data = struct.unpack('{n}h'.format(n=nframes), data)
#    data = np.array(data)
#    return Fs,data,framerate

if __name__ == '__main__':
    Fs = 6000
    T_interval = 1/Fs
    Freq_max = Fs/2
    N=1024
    t = np.arange(0,N-1)*T_interval
    freq = np.linspace(0,Freq_max,N/2+1)
    sig = np.sin(500*np.pi*t)+np.sin(1000*2*np.pi*t)+np.sin(2500*2*np.pi*t)
    freq,amplitude,phase,index,number=freq_amp(sig,N,Fs)
    signature=band_pass(freq,sig,Fs,N)

解决方案

仔细调试下，看看每一步计算和你书上的公式是否一致

解决方案二：

翻书看看你的算法公式是不是对的