Referencia: Ejemplo León-García 4.36 p180
Sea Y=cos(X), y supondremos que X es uniformemente distribuida en el intervalo de (0, 2π]. Y puede ser vista como una muestra de una señal sinusoidal a un instante aleatorio de tiempo que esta uniformemente distribuida en el periodo de la sinusoide.
Para encontrar la pdf, se tiene que la variable -1< y < 1
f_Y(y)=\left.\sum_{k} \frac{f_X(x)}{|dy/dx|} \right|_{x=x_k}La función fX(x) es uniforme, por lo que el área del integral en el intervalo debe ser 1, es decir:
\int_{0}^{2\pi}c dx=1 =(2\pi-0)*c= 2\pi c = 1 c=\frac{1}{2\pi} f_X(x) = \frac{1}{2\pi}De la figura se encuentra que la ecuación tiene dos puntos de intersección o dos soluciones:
x_0 = \cos ^{-1}(y) x_1 = 2\pi- x_0por lo que se evalúa la derivada en cada punto de intersección:
\left.\frac{dy}{dx} \right|_{x_0}= -sin(x_0) = - sin(cos^{-1}(y)) \alpha = \cos^{-1}(x) \sin(\alpha) = y x^2 + y^2 =1 x = \sqrt{1-y^2} \left.\frac{dy}{dx} \right|_{x_0}= -\sqrt{1-y^2} \left.\frac{dy}{dx} \right|_{x_1}= -\sqrt{1-y^2}dado que y0 = y1 = y
Aplicando la ecuación primera:
f_Y(y) = \frac{1}{2\pi \sqrt{1-y^2}} + \frac{1}{2\pi \sqrt{1-y^2}} = \frac{1}{\pi \sqrt{1-y^2}} -1 < y < 1Gráfica del resultado
# pdf de una señal senoidal import matplotlib.pyplot as plt import numpy as np # INGRESO # n=int(input('numero de muestras: ')) n = 200 # PROCEDIMIENTO t = np.linspace(-1,1,n) fy = np.zeros(n,dtype=float) # Evaluar denominador no cero for i in range(0,n,1): denominador = 1-t[i]**2 if (denominador!=0): fy[i] = 1/(np.pi*np.sqrt(1-t[i]**2)) dt = t[1]-t[0] Fy = np.cumsum(fy)*dt # SALIDA plt.plot(t,fy,label='pdf') plt.plot(t,Fy,label='cdf') plt.xlabel('x') plt.legend() plt.show()
Cálculo con un método numérico
Señal Coseno(x), pmf, cdf
Se obtiene barriendo una ventana horizontal estrecha de Δx voltios de ancho, verticalmente a lo largo de las formas de onda y después midiendo la frecuencia relativa de la ocurrencia de voltajes en la ventana Δx.
El eje de tiempo se divide en n intervalos y la forma de onda aparece nΔx veces dentro de estos intervalos en la ventana Δx.
Para observar lo indicado usando python, se presenta el siguiente grupo de instrucciones:
Genera la señal
En un intervalo [a,b] de tiempo, con n muestras en el intervalo y una señal coseno(x) con frecuencia fs Hz.
# Pmf de una senal periódica conocida # definir la señal en la sección senal_t import numpy as np import scipy.signal as signal import scipy.stats as stats import matplotlib.pyplot as plt def senal_t(t, fs): m = len(t) y = np.zeros(m, dtype=float) for i in range(0,m,1): y[i] = np.cos(2 * np.pi * fs * t[i]) return(y) # INGRESO # Frecuencia de la señal en Hz fs = 5 # Rango en tiempo de la señal y muestras a = 0 b = 1 tramos = 2000 # PROCEDIMIENTO muestreo = tramos+1 t = np.linspace(a, b, muestreo) senal = senal_t(t, fs) # SALIDA plt.plot(t,senal) plt.title(' Señal sen(x)') plt.xlabel('t') plt.ylabel('señal') plt.show()
Determinar la función de probabilidad de masa PMF
Se divide el rango de valores posibles de la señal en m intervalos.
Para simplificar el conteo de valores por rango se usa la función scypy.stats.relfreq, que cuenta en m intervalos entre el rango de la señal, dando como resultado las frecuencias relativas a intervalos Δx.
\sum_{i=0}^{N} p_i(x)=1
Si se considera aproximar la función a contínua para observación, será necesario compensar el vector de relativas al dividir por Δx pues el conteo relativo depende del número de divisiones m.
# Pmf de una senal periódica conocida # definir la señal en la sección senal_t import numpy as np import scipy.signal as signal import scipy.stats as stats import matplotlib.pyplot as plt def senal_t(t, fs): m = len(t) y = np.zeros(m, dtype=float) for i in range(0,m,1): y[i] = np.cos(2 * np.pi * fs * t[i]) return(y) # INGRESO # Frecuencia de la señal en Hz fs = 5 # Rango en tiempo de la señal y muestras a = 0 b = 1 tramos = 2000 # PROCEDIMIENTO muestreo = tramos+1 t = np.linspace(a, b, muestreo) senal = senal_t(t, fs) # SALIDA plt.plot(t,senal) plt.title(' Señal sen(x)') plt.xlabel('t') plt.ylabel('señal') plt.show() # Función de Probabilidad de Masa, PMF m = 50 # intervalos en análisis #PROCEDIMIENTO relativa = stats.relfreq(senal, numbins = m ) deltax = relativa.binsize # considerando el deltax, para valor de PMF relcontinua=relativa.frequency/deltax # Eje de frecuencias, por cada deltax senalmin = np.min(senal) senalmax = np.max(senal) senalrango = np.linspace(senalmin,senalmax,m) # SALIDA print('frecuencia relativa:') print(relativa.frequency) print('Rango de Señal') print(senalrango) print('Aproximación a contínua') print(relcontinua) # SALIDA Grafico de PMF plt.subplot(211) plt.bar(senalrango,relativa.frequency, width=deltax*0.8) plt.xlabel('Amplitud de señal') plt.ylabel('PMF') plt.subplot(212) plt.plot(senalrango,relcontinua) plt.ylim(ymin=0,ymax=np.max(relcontinua)) plt.xlabel('Amplitud de señal') plt.ylabel('Aproximación a PDF') plt.show()
frecuencia relativa: [ 0.045 0.034 0.023 0.018 0.017 0.014 0.014 0.012 0.012 0.011 0.011 0.01 0.01 0.01 0.009 0.009 0.009 0.008 0.008 0.009 0.008 0.007 0.008 0.008 0.007 0.008 0.007 0.007 0.007 0.007 ... 0.008 0.007 0.008 0.008 0.008 0.008 0.008 0.009 0.009 0.009 0.009 0.009 0.01 0.01 0.011 0.011 0.012 0.013 0.013 0.015 0.016 0.019 0.023 0.033 0.046] Rango de Señal [-0.99999877 -0.97979676 -0.95959475 -0.93939274 -0.91919073 -0.89898873 -0.87878672 -0.85858471 -0.8383827 -0.8181807 -0.79797869 -0.77777668 -0.75757467 -0.73737266 -0.71717066 -0.69696865 -0.67676664 -0.65656463 ... 0.63636386 0.65656587 0.67676788 0.69696988 0.71717189 0.7373739 0.75757591 0.77777791 0.79797992 0.81818193 0.83838394 0.85858595 0.87878795 0.89898996 0.91919197 0.93939398 0.95959598 0.97979799 1. ] Aproximación a contínua [ 2.22750138 1.68300104 1.1385007 0.89100055 0.84150052 0.69300043 0.69300043 0.59400037 0.59400037 0.54450034 0.54450034 0.49500031 0.49500031 0.49500031 0.44550028 0.44550028 0.44550028 0.39600024 ... 0.39600024 0.44550028 0.44550028 0.44550028 0.44550028 0.44550028 0.49500031 0.49500031 0.54450034 0.54450034 0.59400037 0.6435004 0.6435004 0.74250046 0.79200049 0.94050058 1.1385007 1.63350101 2.27700141]
Determinar la Función de distribucíón acumulada
Considere el caso discreto, solo es necesario ir acumulando los valores de «relativa».
Al usar la aproximación a contínua, la sumatoria debe considerar el área bajo la curva, por lo que los valores se multiplican por Δx antes de acumular los valores de las frecuencias.
# Función distribución acumulada acumulada=np.cumsum(relcontinua*deltax) # Salida CDF # plt.step(senalrango,relcontinua,label='pdf', where='post') plt.step(senalrango,acumulada,label='cdf', where='post') plt.xlabel('Amplitud de señal') plt.title(' Función de distribuión acumulada , CDF') plt.legend() plt.show()
Tarea: realice el ejercicio con otras señales periódicas.
Referencia: Leon W Couch apéndice B p675