Уведомления

Группа в Telegram: @pythonsu

#1 Май 23, 2013 18:14:31

Win95
От:
Зарегистрирован: 2011-11-04
Сообщения: 56
Репутация: +  0  -
Профиль   Отправить e-mail  

Интерполяция, аппроксимация

Доброго времени суток!
Есть некий массив данных по которому строю график.
График получается неравномерным, присутствуют скачки. Возможно ли сгладить график как на рисунке (а лучше свести его к прямой) с помощью numpy или scipy?
Пробовал разобраться в интерполяции:
Scipy
Nympy

Например это:

#a - некий массив
from scipy import interpolate
x = np.arange(0,len(a))
y = a
f = interpolate.interp1d(x, y)
но все время получаю ошибку, хотя длины массивов одинаковые:
    raise ValueError("x and y arrays must be equal in length along "
ValueError: x and y arrays must be equal in length along interpolation axis.

Возможно кто-то сталкивался с похожей проблемой и подскажет куда копать.



Отредактировано Win95 (Май 23, 2013 18:18:06)

Прикреплённый файлы:
attachment xAlLhn.png.pagespeed.ic.9uy0ZlrDiB.png (23,5 KБ)

Офлайн

#2 Май 23, 2013 20:08:42

4kpt
От: Харьков
Зарегистрирован: 2010-11-03
Сообщения: 998
Репутация: +  63  -
Профиль   Отправить e-mail  

Интерполяция, аппроксимация

Приведите весь код. У меня все работает отлчино.

P.S. Если вы хотите построить полином и определить его коэффициенты, то это не то. Вообще не то :)



Офлайн

#3 Май 23, 2013 20:27:09

Win95
От:
Зарегистрирован: 2011-11-04
Сообщения: 56
Репутация: +  0  -
Профиль   Отправить e-mail  

Интерполяция, аппроксимация

Нет проблем. Вот код:
+ вложение это файл необходимый для работы

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy.fft
import pylab
import math
import wave
import numpy as np
import os
import datetime
import scipy
import obspy.signal
from scipy import interpolate
def dependencies_for_myprogram():
    from scipy.sparse.csgraph import _validation
print u'* начинаем перекодирование'
np.set_printoptions(threshold=np.nan)
types = {
    1: np.int8,
    2: np.int16,
    4: np.int32
}
wav = wave.open("record.wav", mode="r")
(nchannels, sampwidth, framerate, nframes, comptype, compname) = wav.getparams()
content = wav.readframes(nframes)
samples = np.fromstring(content, dtype=types[sampwidth])
wav.close()
print u"Определение параметров..."
print u"Ширина выборки", sampwidth
print u"Дискритизация", framerate
print u"Количество частей", nframes
print u"Тип сжатия", comptype
print u"Имя сжатия", compname
print u"Длина файла (секунд)", float(nframes)/framerate
#приведение к нормальному виду матрицы (для записи в файл)
p=samples.reshape(-1, 2)
toshi=len(p)
print u"Количество точек", toshi
#середина массива"
tosh=toshi/2
#отдиление 1 канала от другого (1 канал)
one=np.delete(p, 0, axis=1)
two=np.delete(p, 1, axis=1)
print u"Быстрого преобразование Фурье..."
#быстрое преобразование фурье 
on=(np.fft.rfft(one, axis=0))
tw=(np.fft.rfft(two, axis=0))
#выделяем положительную часть
print u"Выделяем положительную часть..."
#o=((on[0:tosh])*2)
#t=((tw[0:tosh])*2)
o=on
t=tw
#выделяем наибольшие из гармоник
print u"Находим максимальные гармоники..."
oarg=o.argmax(axis=0)
targ=t.argmax(axis=0)
print oarg
print targ
#выделяем максимальные гармоники (сигнал 1)
verh=abs(oarg+100)
niz=abs(oarg-100)
if oarg<100: niz=0
if verh>toshi: verh=toshi    
sig=((o[niz:verh]))
#выделяем 2ю макс гармонику (сигнал 1)
sigl=((o[0:niz]))
try: oa=sigl.argmax(axis=0)
except(ValueError): oa=0
sigld=((o[verh:tosh]))
try:oad=sigld.argmax(axis=0)
except(ValueError): oad=0
loi=oa
loid=oad
if oad>oa: oa=oad
verhl=abs(oa+100)
nizl=abs(oa-100)
if oa<100: nizl=0
if verhl>toshi: verhl=toshi 
print oa
if loid>loi:
    siglo=((sigld[nizl:verhl]))
else: siglo=((sigl[nizl:verhl]))
#выделяем максимальные гармоники (сигнал 2)
up=abs(targ+100)
down=abs(targ-100)
if targ<100: down=0
if up>toshi: up=toshi
signal=((t[down:up]))
#выделяем 2ю макс гармонику (сигнал 2)
signall=((t[up:tosh]))
try: oi=signall.argmax(axis=0)
except(ValueError): oi=0
signalld=((t[0:down]))
try: oid=signalld.argmax(axis=0)
except(ValueError): oid=0
uop=oi
uopd=oid
if oid>oi: oi=oid
upl=abs(oi+100)
downl=abs(oi-100)
if oi<100: downl=0
if upl>toshi: upl=toshi 
print oi
if uopd>uop:
    signallo=((signalld[downl:upl]))
else: signallo=((signall[downl:upl]))
#Обратное преобразование фурье
print u"Обратное проеобразование Фурье..."
q=(np.fft.irfft(sig,axis=0))
w=(np.fft.irfft(siglo,axis=0))
e=(np.fft.irfft(signal,axis=0))
r=(np.fft.irfft(signallo,axis=0))
#Огибающая
a=obspy.signal.filter.envelope(q)
s=obspy.signal.filter.envelope(w)
d=obspy.signal.filter.envelope(e)
f=obspy.signal.filter.envelope(r)
a=(a[10:(len(a)-10)])
s=(s[10:(len(s)-10)])
d=(d[10:(len(d)-10)])
f=(f[10:(len(f)-10)])
x = np.arange(0,len(a))
y = a
f = interpolate.interp1d(x, y)
#Построение графиков
pylab.figure(1)
pylab.subplot(331)
pylab.plot(p)
pylab.xlim (0, 150)
pylab.title (u"Signal")
pylab.subplot(332)
pylab.plot(one, color='green', )
pylab.xlim (0, 150)
pylab.title (u"Signal (1)")
pylab.subplot(333)
pylab.plot(o, color='green',)
pylab.title (u"Range (1)")
pylab.subplot(334)
pylab.plot(sig, color='green')
pylab.title (u"Bandwidth allocation (1.1)")
pylab.subplot(335)
pylab.plot(siglo, color='green')
pylab.title (u"Bandwidth allocation (1.2)")
pylab.subplot(336)
pylab.plot(q, color='green')
pylab.plot(a, color='red')
pylab.title (u"Inverse Fourier transform (1.1)")
pylab.subplot(337)
pylab.plot(w, color='green')
pylab.plot(s, color='red')
pylab.title (u"Inverse Fourier transform (1.2)")
pylab.subplot(338)
pylab.plot(a, color='green')
pylab.title (u"Approximation (1.1)")
pylab.subplot(339)
pylab.plot(s, color='green')
pylab.title (u"Approximation (1.2)")
pylab.figure(1).set_size_inches(24,14)
pylab.savefig('1.png')
pylab.figure(2)
pylab.subplot(331)
pylab.plot(p)
pylab.xlim (0, 150)
pylab.title (u"Signal")
pylab.subplot(332)
pylab.plot(two, color='blue', )
pylab.xlim (0, 150)
pylab.title (u"Signal (2)")
pylab.subplot(333)
pylab.plot(t, color='blue',)
pylab.title (u"Range (2)")
pylab.subplot(334)
pylab.plot(signal, color='blue')
pylab.title (u"Bandwidth allocation (2.1)")
pylab.subplot(335)
pylab.plot(signallo, color='blue')
pylab.title (u"Bandwidth allocation (2.2)")
pylab.subplot(336)
pylab.plot(e, color='blue')
pylab.plot(d, color='red')
pylab.title (u"Inverse Fourier transform (2.1)")
pylab.subplot(337)
pylab.plot(r, color='blue')
pylab.plot(f, color='red')
pylab.title (u"Inverse Fourier transform (2.2)")
pylab.subplot(338)
pylab.plot(d, color='blue')
pylab.title (u"Approximation (2.1)")
pylab.subplot(339)
pylab.plot(f, color='blue')
pylab.title (u"Approximation (2.2)")
pylab.figure(2).set_size_inches(24,14)
pylab.savefig('2.png')
pylab.show()



Отредактировано Win95 (Май 23, 2013 20:28:10)

Прикреплённый файлы:
attachment record.wav (953,2 KБ)

Офлайн

#4 Май 23, 2013 22:25:14

4kpt
От: Харьков
Зарегистрирован: 2010-11-03
Сообщения: 998
Репутация: +  63  -
Профиль   Отправить e-mail  

Интерполяция, аппроксимация

Что Вы в итоге хотите получить?
- Красивый график.
- Полином необходимой степени.
Так будет легче пилить Ваш файл…

P.S. Monkey?



Отредактировано 4kpt (Май 23, 2013 22:26:01)

Офлайн

#5 Май 23, 2013 22:37:24

Win95
От:
Зарегистрирован: 2011-11-04
Сообщения: 56
Репутация: +  0  -
Профиль   Отправить e-mail  

Интерполяция, аппроксимация

я хочу получить верный график. мне необходимо произвести действие над огибающей

a=obspy.signal.filter.envelope(q)
s=obspy.signal.filter.envelope(w)
d=obspy.signal.filter.envelope(e)
f=obspy.signal.filter.envelope(r)
(аппроксимировать, интерполировать, сделать что-то еще) и в конечном итоге привести ее вид к прямой с которой я буду работать дальше



Офлайн

#6 Май 24, 2013 01:06:07

4kpt
От: Харьков
Зарегистрирован: 2010-11-03
Сообщения: 998
Репутация: +  63  -
Профиль   Отправить e-mail  

Интерполяция, аппроксимация

Ох и код… Мда…
Запустил и получил результат. Это не то, что Вы ожидали.
Если я правильно понял, то нужно найти для a, s, d, f - полиномы первой степени, т.е. коэффициенты уравнения y = ax + b (найти коэффициенты a, b). При получении коэффициентов необходимо минимизировать ошибку межту результатами уравнения в точках и эмпирическими данными.

Я все правильно понял? Если правильно, то завтра напишу…

Единственный допвопрос: что является для этих данных осью х - просто ряд от 0 до их размера?



Офлайн

#7 Май 24, 2013 07:46:05

Win95
От:
Зарегистрирован: 2011-11-04
Сообщения: 56
Репутация: +  0  -
Профиль   Отправить e-mail  

Интерполяция, аппроксимация

Вы верно поняли. Изначально осью х является время (это же wav файл), т.е от точки до точки проходит 0,05 секунд (зависит от дискретизации), однако оси времени (и амплитуды) я не переводил, по этому Вы правы - осью х является длина выбранного отрезка.



Офлайн

#8 Май 26, 2013 00:49:58

4kpt
От: Харьков
Зарегистрирован: 2010-11-03
Сообщения: 998
Репутация: +  63  -
Профиль   Отправить e-mail  

Интерполяция, аппроксимация

Код еще нужен? Я просто выпал на пару дней…



Офлайн

#9 Май 26, 2013 16:17:11

Win95
От:
Зарегистрирован: 2011-11-04
Сообщения: 56
Репутация: +  0  -
Профиль   Отправить e-mail  

Интерполяция, аппроксимация

4kpt
Код еще нужен? Я просто выпал на пару дней…
Да, еще нужен. Пытался сам разобраться с полиномами , нашел это.
Сейчас разбираюсь.



Офлайн

#10 Май 26, 2013 16:42:38

4kpt
От: Харьков
Зарегистрирован: 2010-11-03
Сообщения: 998
Репутация: +  63  -
Профиль   Отправить e-mail  

Интерполяция, аппроксимация

Проблема метода заключается в определении этих весовых коэффициентов. Я предпочитаю классический способ апроксимации…



Офлайн

Board footer

Модераторировать

Powered by DjangoBB

Lo-Fi Version