Android-приложение для поиска дешевых авиабилетов: play.google.com
Главная -> Справочник по алгоритмам

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 [40] 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78

ввод всех отсчетов г/,-. Если память ПЭВМ достаточна, можно предусмотреть хранение yi с помощью массива Y (/), содержащего действующие (т. е. У1ф0) отсчеты. Тогда, предусмотрев в п. 3 алгоритма вместо ручного ввода yi циклический вызов отсчетов из массива Y (/), получим вариант комбинированного или последовательно-параллельного спектрального анализа. В приведенной ниже программе 5.2 реализован этот метод. Кроме того, в ней предусмотрен ввод опорной частоты f\ = Fl. Это частота первой гармоники при анализе периодических функций или частота, задающая масштаб сетки частот kf\ при анализе финитных функций (т. е. в последнем случае им придается характер периодических функций с частотой повторения f,). Коррекция Mk и фазового сдвига в этой программе не предусмотрена (но может быть легко введена). Программа 5.2.

Значения М4 и S4 недостоверны ввиду их малости (точные значения Mi==0 и 54 = 0). Поэтому случайные значения Q4 из выдачи исключаются. Для t - 0,5-l0~ по программе 5.2 получим 1/(0 = 1,07593, а для / = = 3-10-" получим у (0=0,00852.

Спектральный анализ функций с ограниченным спектром выполняется также непосредственно с помощью дискретного преобразования Фурье над комплексными числами. Так, если задан массив и,-=А-, + уТ, ( = 0, 1, 2, ..., /V-1), то прямое дискретное преобразование Фурье дает массив комплексных чисел

.V-I 1 = 0

-i {2nki/N)

Соответственно массив чисел и* при обратном дискретном преобразовании Фурье обра-

10 PRINTПОСЛЕДОВАТЕЛЬНО-ПАРАЛЛЕЛЬНЫЙ СПЕКТРАЛЬНЫЙ АНАЛИЗ

15 PRINTH ТРИГОНОМЕТРИЧЕСКАЯ ИНТЕРПОЛЯЦИЯ-ЭКСТРАПОЛЯЦИЯ

20 ШРиТЗАДАйТЕ ОБЩЕЕ ЧИСЛО ОТСЧЕТОВ ФУНКЦИИ V<T) N=N

30 INPUTЗАДАЙТЕ НОМЕРА НАЧАЛЬНОГО И КОНЕЧНОГО ОТСЧЕТОВ IS,IF U/G

40 LETG=6-U!DIM V<:G>/A<N2)/B<:Nx2)!LETR=2iif#PIN

50 FOR 1=0 TO G:PRINT!3.0!ВВЕДИТЕ Vl+U:INPUT VtOsHEXT I

60 INPUTЗАДАЙТЕ ЧИСЛО СОСТАВЛЯ»ЩИХ СПЕКТРА М=М

70 INPUTЗАДАЙТЕ НОМЕР ПЕРВОЙ СОСТАВЛЯЮЩЕЙ СПЕКТРА KS=S

75 INPUTЗАДАЙТЕ ЧАСТОТУ F1=F

80 LETA=0sLETB=0!LETC=0!FOR К=1 ТО М

90 LETA<:K)=0!LETB<K)=0:NEXT к

100 FOR 1=8 ТО G!LETA=A+V<;n:LET2=R«<I+U)

110 FOR К=1 ТО М:ЬЕТи=2ж<К+8-1)

128 LETA<:K)=A<K>+V<I)«COS<W>!LETB<K)=B<K>-bV<I>*SINaJ)

130 NEXT к:NEXT l!LETM0=AN:LETS0=MeF

140 PR INT! F 1.9!ПОСТОЯННАЯ СОСТАВ ЛЯК1ЩАЯ М0=М0

158 PRINTСПЕКТРАЛЬНАЯ ПЛОТНОСТЬ НА НУЛЕВОЙ ЧАСТОТЕ 88=88

160 FOR К=1 ТО M!LETU=SGR<A<K>«A<K>+BCK>«B<K>)

178 LETS=-ACS<A<K)U):IF БСКХ0 THEN LETS=-G

180 LETO=K+S-l!LETB=B+U«U:IF K+S=2 THEN LETC=U

190 PRINT!3.0!АМПЛИТУДА M0!F1.9!=U«2N

200 PRIHT!3.0!СПЕКТРАЛЬНАЯ ПЛОТНОСТЬ S0!F1.9! = U.N.F

210 PR INT! 3.0!ФАЗОВЫЙ СДВИГ e0!F1.9!=DEG<:e>

220 NEXT K:IFS<>1 THEN 60

230 PR I NTКОЭФФИЦИЕНТ ГАРМОНИК КГ=!F1.9!SQR<:B-CiifC>C 240 INPUT1-БУДЕТ ВЫЧИСЛЯТЬСЯ V<T>,e-HE БУДЕТ ?D 250 IF DOl THEN 60

260 ШРиТЗАДАйТЕ ВРЕМЯ T=T!LETV=e!LETP=2*#PI*TiifF 270 FOR K=l TO M!LETW=P*KsLETV=V+A<:K)iifCOS<W>+B<K>«SIN<:W> 288 NEXT K:LETV=M0+Viif2/N:PRINTЗНАЧЕНИЕ V<T)=!F1.9!V 290 GOTO 260:END

Пример. Для примера к программе 5.1 (Л/ = 32, (/„-(/7 = 1, =250ООО), задав М = 5 и УИ5=1, получим (даны 5 цифр после запятой):

.МО = 0,25

Ml =0,45088

.М2 = 0,32036

.МЗ = 0,15224

М4 = 2,23345-10-*

М5 = 0,09375

Кг = 0,81369

Q1 = -39,375 Q2=-78,75 Q3= -118,125 Q5=-16,875

SO=MO- S1=9,D176-10-

52 = 6,40729-10-

53 = 3,0449-10-

54 = 4.466892-10- 55=1,87503-10-

зует массив комплексных чисел

M/N)

* = 0

Если Ui=Y (ti), причем Re u, = X,= У (?,). a lmu,= y, = 0, TO прямое .преобразование Фурье дает коэффициенты С(,=Ли) -уВи) = = vi,/N усеченного тригонометрического ряда Фурье в комплексной форме:

y(t)=l iai."<*+"*.

где \Ск\=л1х1-{-П, ф*=-агссо5 (Х,/С*1), если У/,>0, и ф(,= -фс, если У*<0. Пер-



вые Hjl Значений Сь соответствуют положительным частотам- \ = Щ\, остальные - отрицательным.

Соответственно, Hjvien частотный спектр в виде Л значений СсХ + уТ* и подвергнув его обратному преобразованию Фурье, получим Л/ значений </, (ti)/N. Из них первые Л 2 значений соответствуют области времен /0, остальные - /<0.

Указанные дискретные преобразования вытекают из прямого преобразования Фурье

S (ш)= \ y(t) €-•dt

- оо

и ИЗ обратного преобразования Фурье

выполняемых методиками приближенного численного интегрирования для функций у (t), определенных на конечном промежутке вре-

точно вычислять 2М раз, поскольку они периодически повторяются. Кроме того, вычисления можно проводить по рекуррентным формулам (см. подробно в [5]). Выделив четную и нечетную части разложений, удается уменьшить время вычислений еще вдвое; в результате БПФ требует проведения (Л 2) log2 N операций вместо Л/* (при Л/=1024 это уменьшает число операций более чем в 200 раз!). Для ПЭВМ характерны значения Л/64-256, а для профессиональных ПЭВМ - в несколько раз больше.

Выбор N из ряда 2 (2, 4, 8, 16 и т. д.) является определенным неудобством БПФ. Однако оно компенсируется заметным уже при Л/= 8 или N=16 сокращением времени вычислений. Программы БПФ обычно позволяют проводить как прямое, так и обратное дискретное преобразование Фурье над массивами ui=A;-f-/T; и vii = Xi,-\-iY;, обеспечивая преобразование у,-» Yn, и наоборот.

Программа 5.3.

10 prihtпрямое и обратное быстрое пре0бра30б**!ие фурье

£0 inputдля n=2M введите M=M:letn=int<:2M+. OsDIM хсн-п, VCN-O

30 inputзадайте -1 при прямом бп* и 1 при обратном d

35 inputвведите номера начального и конечного отсчетов IS,If 10,11

40 for 1=0 то 11-10!рр1нт!3.0!для i=I+I0

50 inputвведите XbVI XCDVvCDsNEXT I

60 for l=l to M!lete=intc£"<M+l-l>+.l)!letf=e/2!letu=l!letu=0

70 let2=#pif:letc=c0s<:2)!lets=d!iiSIN<:2>!f0r j=l to f

80 for i=j to n step e:leto=I+f-l:letp=x<i-l>+X<0)!letg=Va-l>+v<0>

90LETR=x<:i-l>-x<o)!lett=v<:i-l>-v<o>:letx<:o>=r*u-TiifU

100 letv<:0>=TiifU+r*u!letx<I-l>=p!letv<I-l>=g!next i

110 ЕЕТЫ=ижС-иж8:иЕТи=ижС+ижЗ!ЕЕТи=Ы:нехт JsNEXT l

120 letj=l:for 1=1 to n-l:if i>=j then 150

130 letj1=j-1!leti 1=1-1:letp=x<:л >sletg=v<Jl):letx<Jl)=X<11)

140 LETV<:j1>=V<I1)!LETX<:I1)=P:LETV<I1)=Q

150 letk=n/2

160 if k>=j then 180

170 letj=j-K:letk=k2!60t0 160

180 LETJ=J+K:NEXT I

190 if d=-l then printрезультаты прямого бп*!бото 240 200 printрезультаты обратного бпф:for к=0 to n-l 210 letx<k>=x<k>N!Letv<:k)=v<:k)N

220 PRINT!3.0!k=Ki!f1.5! x=x<k>i v=v<k) 230 next к:stop"

240 for k=0 TO n-l!leta=ssr<x<:k)«x<k)+v<:k)iifV(K)) .

250 letq=0:if a=0 тжн 270

260 letq=acs(x<k>a>:if v<k><0 then LETe=-e

270 print !3.e!k=k;!F1.5! х=х<:к> v=v<k);

280 print M=a«2/-Ni g=de6<s)!next KsEND

мени, и зависимостей S (w) .с ограниченным спектром.

Дискретное преобразование Фрье требует числа операций порядка Л/" и ведет к большим затратам машинного времени уже при N> I0-20. Кроме того (при one- . рациях с комплексными числами), приходится 2N раз вычислять тригонометрические функции (sin н cos) - операции по разложениям их в ряд выполняются на ЭВМ довольно медленно.

Быстрое преобразование Фурье (БПФ) позволяет резко уменьшить время проведения прямого и обратного дискретных преобразований Фурье. При БПФ N выбирается из

условия N = 2", где /М = 1, 2, 3.....При этом

значения тригонометрических функций доста-

Пример. В программе 5.3 предусмотрен ввод только ненулевых отсчетов. Например, для вычисления спектра прямоугольного импульса, заданного при Л/= 32 (М = 5) восемью единичными отсчетами Y (t) Хо, Х,, Xi, ..., и остальными нулевыми (все у; = 0), задаем номер начального отсчета /5=0 и конечного /£=7 и указываем тип преобразования. При прямом дискретном преобразовании Фурье получаем спектральный состав (даны первые 6 значений из таблицы результатов, выводимой на экран дисплея):

К=0 Х=8 Y=0

К=1 Х = 5,57659 Y=-4,57659

К = 2 Х=1 ¥=-5,02734



К=3 Х=-1,14828 Y=-2,1482e

К=4 Х=0 Y=0

К=»=5-Х== 1,43543 Y=-1.43543

М==0,25 Q = 0

М=0,45088 Q = - 39,75

М=0,32036 Q=-78,75

М = 0,15224 Q=-118,125

М=0 Q=0

М = 0,09375 Q=-16,875

Если дать команды D = 1 GOTO 60, то после прямого будет проведено обратное дискретное преобразование Фурье. Без учета малых погрешностей округления, получим в итоге значения Хс~-Хт1. Xs-Xi-O и УоУз! = 0 (т. е. исходный массив ui). Время преобразования, при М = 5 около 2,5 мин.

Повышение точности БПФ достигается специальной обработкой входных или выходных данных. Обработка входных данных заключается в умножении отсчетов на весовые коэффициенты Wi, подобранные так, чтобы обеспечить заданную аппроксимацию Y (t) в ходе интегрирования. Значения Wi~\ соответствуют ступенчатой аппроксимации У (t) и интегрированию методом прямоугольников. Если 1,= 1/2, 1, I, ...

1, 1/2, то У (О в ходе БПФ аппроксимируется кусочно-линейной функцией, т. е. интегрирование происходит методом трапеций. При 1Г, = 1/3, 4/3, 2/3, 4/3, 2/3, ...

4/3, 1/3 аппроксимация Y (t) будет параболической, т. е. интегрирование происходит по методу Симпсона. При более сложных последовательностях W, возможна аппроксимация У (i) полиномами более высоких степеней.

Применим этот метод к примеру для программы 5.3. Задад11м 9 отсчетов Х()=1/2, Х,Хг=1 и А:8=1/2 (все Y„-Yi=0). С помощью программы 5.3 получим (первые 6 строк таблицы выходных данных):

К = 0 Х=8 Y = 0

К = 1 X = 5,07659 Y = - 5,07659

К = 2 Х=-7-10- Y=-5.02734

К = 3 1,64828 Y=-1,64828

К -4 Х=0 Y=0

К = 5 Х = 0,93543 Y=-0,93543

М=0,25 Q = 0

М = 0,44871 Q=-45

М = 0,31421 Q=-90

М=0,14569 Q= -135

М=0 Q=0

М = 0,09268 Q-45

При параболической аппроксимации У (/) задаем ее значениями Хо-1/3, Ai=4/3, Лг = 2/3, 5=4/3, А, = 2/3, As = 4/3, А-с = = 2/3, А7 = 4/3, х8=1/3. х9-4-Аэ1==0 и Уо-;-Хз1=0. Тогда получим:

К=0 Х=8 Y=0

К=1 Х=5,0930 Y=-5,09.30

К = 2 Х=6.6-10- Y=-5,09364

К = 3 Х=-1,69884 Y=-l,69884

к=4 Х=0 Y=0

К = 5 Х= 1,02452 Y= -1,02452

М=.0 Q=.0

.М = 0,45016 Q=-45

М = 0,31835 Q=-90

М = 0,15016 Q==-135 •

М=0 Q=0 • .

М=0,0906 Q=-45

Для данного примера точные результаты вычисляются по формулам

X* = sin {nk/2) N/2nk.

2t=-(l-cos(nfe/2)) N/2nk,

M*=2Vr+y!/N.

Q*=-arccos (Xt/VF+yf)

И. равны:

K=0 X=8 Y=0

K=l X = 5,09296 Y=-5,09296

K = 2 X = 0 Y=-5,09296

K = 3 X=-1,69765 Y = -1,69765

K=4 X=0 Y=0

K = 5 X= 1,01859 Y=-1,01858

M = 0,25 Q = 0

М=0,450Г6 Q=-45

M = 0,31831 Q=-90

M = 0,15005 Q= -135

M=0 Q=0

M = 0,09003 Q=-45

Сравнение приведенных данных наглядно иллюстрирует повышение точности БПФ по мере увеличения порядка аппроксимирующего у (О полинома. Казалось бы. что для прямоугольного импульса ступенчатая аппроксимация будет наилучшей. Однако при БПФ производится приближенное численное интегрирование не самой зависимости у {t), а ее произведения на быстроосииллирующие множители. Поэтому результирующая погрешность интегрирования определяется не только погрешностью аппроксимации у (t), ио и всей подынтегральной функции. В следующем разделе описан метод спектрального анализа, в котором по точным формулам интегрируется произведение у (t) на осциллирующие множители (при этом простейшая кусочно-линейная аппроксимация у {t) с разрывами позволяет получить точные результаты спектрального анализа).

Повышение точности БПФ путем обработки выходных данных сводится к умножению Xk и Yk на множитель Ki.= = sm (nk/N)/(nk/N) и Введению поправки на фазовый сдвиг Д(()(, = лА/Л/. В этом случае для рассмотренного примера будем иметь точные значения Хь, Yk. Мк и qjt. Разумеется, нельзя использовать одновременно оба способа повышения точности БПФ, поскольку оии являются альтернативными вариантами одного и того же подхода - улучшение аппроксимации У (t).

Эффект Гиббса наблюдается при тригонометрической интерполяции усеченным рядом Фурье функций с разрывами (в частности, прямоугольного импульса). Он заключается в осцилляциях расчетной зависимости ряда Фурье, амплитуда которых может



0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 [40] 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78



0.0369
Яндекс.Метрика