|
Главная -> Справочник по алгоритмам 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.011 |
|