Форум сайта python.su
0
Доброго времени суток!
Есть некий массив данных по которому строю график.
График получается неравномерным, присутствуют скачки. Возможно ли сгладить график как на рисунке (а лучше свести его к прямой) с помощью 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)
Прикреплённый файлы:
xAlLhn.png.pagespeed.ic.9uy0ZlrDiB.png (23,5 KБ)
Офлайн
63
Приведите весь код. У меня все работает отлчино.
P.S. Если вы хотите построить полином и определить его коэффициенты, то это не то. Вообще не то :)
Офлайн
0
Нет проблем. Вот код:
+ вложение это файл необходимый для работы
#!/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)
Прикреплённый файлы:
record.wav (953,2 KБ)
Офлайн
63
Что Вы в итоге хотите получить?
- Красивый график.
- Полином необходимой степени.
Так будет легче пилить Ваш файл…
P.S. Monkey?
Отредактировано 4kpt (Май 23, 2013 22:26:01)
Офлайн
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)
Офлайн
63
Ох и код… Мда…
Запустил и получил результат. Это не то, что Вы ожидали.
Если я правильно понял, то нужно найти для a, s, d, f - полиномы первой степени, т.е. коэффициенты уравнения y = ax + b (найти коэффициенты a, b). При получении коэффициентов необходимо минимизировать ошибку межту результатами уравнения в точках и эмпирическими данными.
Я все правильно понял? Если правильно, то завтра напишу…
Единственный допвопрос: что является для этих данных осью х - просто ряд от 0 до их размера?
Офлайн
0
Вы верно поняли. Изначально осью х является время (это же wav файл), т.е от точки до точки проходит 0,05 секунд (зависит от дискретизации), однако оси времени (и амплитуды) я не переводил, по этому Вы правы - осью х является длина выбранного отрезка.
Офлайн
63
Код еще нужен? Я просто выпал на пару дней…
Офлайн
0
4kptДа, еще нужен. Пытался сам разобраться с полиномами , нашел это.
Код еще нужен? Я просто выпал на пару дней…
Офлайн
63
Проблема метода заключается в определении этих весовых коэффициентов. Я предпочитаю классический способ апроксимации…
Офлайн