Aaaaaallora.... devo aver fatto un po' di confusione mischiando analogico e digitale.

Il codice corretto dovrebbe essere questo (devo solo studiare perchè con il periodogram devo usare 10 anzichè 20):

from __future__ import unicode_literals


from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
import datetime
from obspy.core import read
import obspy
from obspy import *
from obspy.core import utcdatetime
#from obspy.signal import bandpass
#from numpy import mean, sqrt, square, arange
import operator
import sys
from obspy.signal.trigger import classic_sta_lta
from obspy.signal.trigger import carl_sta_trig
from scipy.signal import bessel, lsim
import pylab
import copy
from scipy import signal
import os
import time
import math
from scipy import integrate
from obspy import read
#import plotly.plotly as py
#import plotly.graph_objs as go


from obspy.signal.trigger import classic_sta_lta, plot_trigger
#intanto vado a leggere la traccia in SAC
from math import*



# creating I order, 10 Hz, low pass Butterworth filter
fc = 10 # Hz
wc = 2*np.pi*fc # rad/s
b, a = signal.butter(1, wc, 'low', analog=True)    # mi ritorna un sistema Num(s)/Den(s)

# verifica della funzione di trasferimento
w, h = signal.freqs(b, a)	# ideal frequency response of filter
w = 0.5*w/np.pi				# rad/s to Hz
amp=20*np.log10(np.abs(h))	# amplitude to dB

plt.plot(w, amp)
plt.xscale('log')
plt.xlabel('Frequency Hz')
plt.ylabel('Attenuation dB')
plt.grid()
#plt.ylim(-60, 0)
plt.show()


# creazione della base dei tempi
time=80.0	# durata del segnales
Fs=1000	# frequenza di campionamento Hz
t = np.linspace(0, time, (time*Fs)+1)

# let's make some noise
mu, sigma = 0, 0.1
noise = np.random.normal(mu, sigma, len(t))  # noise

# simulazione della risposta utilizzando lsim
u = noise
# U e T sono i vettori dei valori e la scala dei tempi
tout, yout, xout = lsim((b, a), U=u, T=t)

# passaggio al dominio della frequenza
f, psd = signal.periodogram(u, Fs)
psd_noise=psd[1:]

f, psd = signal.periodogram(yout, Fs)
psd_signal=psd[1:]

# removing first element due related to 0 Hz (for log scale)
f=f[1:]


# ratio filtered/unfiltered signal
legenda=[]
tf=psd_signal/psd_noise
#plt.plot(f,20*np.log10(tf))
plt.plot(f,10*np.log10(tf))

legenda.append("real")
plt.plot(w, amp)
legenda.append("ideal")
plt.xscale('log')
plt.xlabel('Frequency Hz')
plt.ylabel('Attenuation dB')
plt.legend(legenda)
plt.grid()
#plt.ylim(-60, 0)
plt.show()

A meno di errori di copia e incolla dovrebbe funzionare....

Buongiorno a tutti, sono un novizio di Python e vi chiedo cortesemente aiuto per questo quesito:


Voglio verificare la risposta in frequenza di un filtro quindi prendo un rumore bianco e lo filtro per poi fare il rapporto spettrale e sovrapporre la risposta in frequenza nominale del filtro.


Ovviamente il risultato non coincide e non riesco a capire dove sbaglio....




1) creazione del filtro passa-basso Butterworth I ordine, 10 Hz

# creating I order, 10 Hz, low pass Butterworth filter
fc = 10 # Hz
wc = 2*np.pi*fc # rad/s
b, a = signal.butter(1, wc, 'low', analog=True)    # mi ritorna un sistema Num(s)/Den(s)
num=b
den=a

w, h = signal.freqs(b, a)    # ideal frequency response of filter
w = 0.5*w/np.pi                # rad/s to Hz
amp=20*np.log10(np.abs(h))    # amplitude to dB

'''
# ideal filter frequency response dB vs Hz
plt.plot(w, amp)
plt.xscale('log')
plt.grid()
plt.xlabel("Frequency Hz")
plt.ylabel("Amplitude dB")
plt.show()
'''

2) ricavo la risposta impulsiva del filtro da usare con la convoluzione nel dominio del tempo per filtrare un generico segnale

system=signal.lti(num, den)
t, impulse_resp = signal.impulse(system)
impulse_resp = impulse_resp[1:]     # removing first element due it's always 0.0
t = t[1:]

print("Risposta impulsiva:")        # filter impulse response (kernel of filter in time domain)
print(impulse_resp)

# normalizing kernel to 1.0
kernel=len(impulse_resp)
impulse_resp=impulse_resp /sum(impulse_resp)    # vado a normalizzare il kernel per il filtro
print(sum(impulse_resp))

# now I can use impulse response as kernel in time domain

3) creo un rumore bianco (ricavo lo spettro) e lo filtro con la convoluzione nel dominio del tempo (e rivavo lo spettro anche di questo)

# let's make some noise
mu, sigma = 0, 0.1
sample = 100000            # 100k samples
Fs=1000                    # SPS sampling frequency
x = np.arange(sample)
noise = np.random.normal(mu, sigma, len(x))  # noise

# time domain to frequency domain (noise)
#f, psd = signal.periodogram(noise, Fs)
f, psd = signal.welch(noise, Fs, scaling='spectrum')

# removing first element due related to 0 Hz (for log scale)
f=f[1:]
psd=psd[1:]


# filtering in time domain using convolution and filter impulse response
filtrato=np.convolve(noise, impulse_resp)    # filtering noise
filtrato=filtrato[:len(filtrato)-len(impulse_resp)+1]    # equalizing length due to kernel dimension

# time domain to frequency domain for filtered noise
#f, psd2 = signal.periodogram(filtrato, Fs)
f, psd2 = signal.welch(filtrato, Fs, scaling='spectrum')

# as above
f=f[1:]
psd2=psd2[1:]

4) rapporto degli spettri, plot della risposta in frequenza e confronto con la risposta in frequenza nominale del filtro:

# ratio filtered/unfiltered signal
legenda=[]
tf=psd2/psd
plt.plot(f,20*np.log10(tf))
legenda.append("real")
plt.plot(w, amp)
legenda.append("ideal")
plt.xscale('log')
plt.xlabel('Frequency Hz')
plt.ylabel('Attenuation dB')
plt.legend(legenda)
plt.show()

Magari voi che siete più esperti riuscite a capirlo subito...
Grazie in anticipo.