Перейти из форума на сайт.

НовостиФайловые архивы
ПоискАктивные темыТоп лист
ПравилаКто в on-line?
Вход Забыли пароль? Первый раз на этом сайте? Регистрация
Компьютерный форум Ru.Board » Компьютеры » Прикладное программирование » Закладки » Вопросы программирования на FORTRAN (ФОРТРАН)

Модерирует : ShIvADeSt

 Версия для печати • ПодписатьсяДобавить в закладки
Страницы: 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 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330

Открыть новую тему     Написать ответ в эту тему

akaGM

Platinum Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
Обсуждаются все вопросы, связанные с программированием на ФОРТРАН, как общего так и конкретного характера.
Постарайтесь дать как можно больше информации о возникшей проблеме -- это в конце концов в ваших же интересах чтобы вам помогли...

прежде чем просить помощи в задании
платное решение задач

ресурсы этого топика
ссылка на подборку ресурсов, собранных посетителями этого форума
 


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

Всего записей: 24766 | Зарегистр. 06-12-2002 | Отправлено: 18:11 14-01-2007 | Исправлено: akaGM, 15:26 15-05-2024
akaGM

Platinum Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
XPEHOMETP
"не нужна" -- не значит, что я стремлюсь к
не впадай в крайности...
 
мне надо:
1) любая длина входных данных
2) простота реализации
всё остальное вторично

Всего записей: 24766 | Зарегистр. 06-12-2002 | Отправлено: 15:22 09-04-2011
akaGM

Platinum Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
вопрос знатокам Фурье-преобразования
 
почему во всех реализациях оно считается только на входной сетке, типа fft(src, dest, sizeN),
или здесь речь идёт только о коэффициентах?
 
в общем кто виноват и что делать...
 
-----
это не оффтоп, если посмотреть на подсабж топа и Fourier очень похоже по написанию на Fortran :)

Всего записей: 24766 | Зарегистр. 06-12-2002 | Отправлено: 12:04 11-04-2011 | Исправлено: akaGM, 12:06 11-04-2011
Vskazka

Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
akaGM

Цитата:
почему во всех реализациях оно считается только на входной сетке, типа fft(src, dest, sizeN),
или здесь речь идёт только о коэффициентах?
 
в общем кто виноват и что делать...  

 
Если, я правильно понял вопрос (он несколько экстравагантно задан) то очень коротко так
 
Вообще-то, fft - это не преобразование Фурье, в собственном смысле, Там оно континуальное. В теховской нотации это \int_{-\infty}^{\infty}f(x)\exp(-i \xi x) dx.  FFT имеет дело с  конечным рядом Фурье (конечная сумма) вида \sum_{n=-N}^N a_n \exp(i*n*k*s). "Быстрое преобразование Фурье" fft преобразует один, вообще говоря комплексный, вектор на другой вектор такой-же длины. По определенным правилам.  То есть решается задача:  если нам дан вектор b_k, то как его представить в виде  b_k=\sum_{n=-N}^N a_n exp(i*n*k*s). Другое дело, что из FFT можно слепить бывает (и часто лепят) и вычисление собственно преобразования Фурье.
Но это отдельная песня.  
   

Всего записей: 382 | Зарегистр. 24-11-2003 | Отправлено: 12:48 11-04-2011 | Исправлено: Vskazka, 12:50 11-04-2011
XPEHOMETP

Silver Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
akaGM
A в чем проблемы? "на входной сетке" - в смысле, результаты разложения выводятся во входном массиве, который затирается? Дык, экономия памяти, алгоритмы-то эти еще при царе Горохе разработаны! Не хотите затереть исходные данные - делайте их копию и посылайте ее в разложение по Фурье. Или используйте упомянутые мной программы ezfftf + ezfftb, у них сделан вывод коэффициентов в отдельные массивы.
 
Добавлено:
Vskazka
Смыслов понимания преобразования Фурье - целая куча, по крайней мере три. Конечно, тут идет речь о дискретном преобразовании Фурье, но это уже оффтоп.

Всего записей: 2485 | Зарегистр. 21-06-2005 | Отправлено: 12:54 11-04-2011
akaGM

Platinum Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
да я просто не понимаю, почему
 
fft(вход, выход, n)
 
sizeof(выход)=sizeof(вход)=n
 
 
мне надо получать преобразование того размера, кот. я бы сам задавал...
короче
 
до сих пор я считал Ф простым интегралом в лоб, имея дискретные данные F(xi) размера N
и получал массив Ф(zj) размера M
 
Ф(z) = Integral { F(x)*cos(x*z) dx }
 
разве не оправдана такая задача?
имеем сигнал в виде N-временных отчётов,
хотим получить спектр из М-"частотных" отсчётов
 
---
в общем зря я плохо учился...

Всего записей: 24766 | Зарегистр. 06-12-2002 | Отправлено: 13:14 11-04-2011
Vskazka

Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
XPEHOMETP
Преобразование Фурье (просто, как таковое) - одно.  Дискретное преобразование Фурье - это дискретное преобразование Фурье. Быстрое преобразование Фурье  - это частный случай дискретного. Собственно, черный чай с сахаром и лимоном - это именно черный чай с сахаром и лимоном (самостоятельный напиток), но не черный чай.  
Согласен это офтоп, но из за недоговоренности о терминах масса казусов случается.  
Впрочем, я, похоже не очень правильно и вопрос понял...    
 
 
Добавлено:
akaGM

Цитата:
до сих пор я считал Ф простым интегралом в лоб, имея дискретные данные F(xi) размера N
и получал массив Ф(zj) размера M
 
Ф(z) = Integral { F(x)*cos(x*z) dx }  

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

Всего записей: 382 | Зарегистр. 24-11-2003 | Отправлено: 13:17 11-04-2011 | Исправлено: Vskazka, 13:26 11-04-2011
XPEHOMETP

Silver Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
akaGM

Цитата:
мне надо получать преобразование того размера, кот. я бы сам задавал...

Из штанов не выпрыгнешь... В сторону роста частот - там теорема Найквиста задает предел. Все что свыше - то от лукавого. В сторону уменьшения - так при желании потом можно отбросить все ненужное, что программа выдаст.
 
Кроме того, у Вас разложение только по косинусам, это далеко не во всех случаях прокатит...

Всего записей: 2485 | Зарегистр. 21-06-2005 | Отправлено: 13:25 11-04-2011 | Исправлено: XPEHOMETP, 13:27 11-04-2011
akaGM

Platinum Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
хорошо, спрашиваю в последний раз... и больше спрашивать не буду :)
 
могу ли я (как?), имея набор n дискретных данных Фi в диапазоне T1 .. Tn
при помощи DFT посчитать _наперёд заданное количество_ m Fj точек в диапазоне W1 .. Wm
 
можно на основе того же горячо любимого здесь FFTPack
только пожалуйста на пальцах, без Найквистов-Котельниковых-Шеннонов, частот среза и нижне-верхних фильтров, по рабоче-крестьянски...

Всего записей: 24766 | Зарегистр. 06-12-2002 | Отправлено: 14:06 11-04-2011 | Исправлено: akaGM, 14:08 11-04-2011
XPEHOMETP

Silver Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
akaGM
Гым! Дискретное преобразование Фурье имеет свои реальные пределы. Это примерно как в случае решения системы линейных уравнений. Фактически, имеем N исходных точек, по ним можем определить N значений переменных. И ни на одну больше. Тут - то же самое. Фактически, мы записываем некие данные в виде многочлена, составленного из ортогональных функций - синусов-косинусов от разных частот. Имеем N точек - получаем N коэффициентов многочлена, и ни на одну больше. Половина из них - коэффициенты при синусах, другая половина - при косинусах. И того - частот получится N/2, что как раз отвечает теореме Найквиста. Что логично.
 
Да, тут есть тонкости. Кроме синусов-косинусов в разложении Фурье есть еще свободный член, среднее арифметическое из всех данных. Но на пределе Найквиста, при N/2 (задавая волну синуса-косинуса всего двумя точками), мы на самом деле не сможем одинаково хорошо определить и синус, и косинус. Задать волну косинуса двумя точками (в начале и в середине цикла) вообще не реально - это будут нули (по определению косинуса), и амплитуда волны тогда совсем не определена. Значит, коэффициент при косинусе на частоте N/2 по определению - ноль, и искать надо только N переменных из N точек данных. Это радует.
 
Тонкостей, на самом деле, еще больше, ибо выше полагалось, что число N - четное. Для нечетного N все немного сложнее, но не будем туда влезать.
 
Короче. Ровно так же, как решая N линейных уравнений относительно N переменных, можно получить ровно N их (переменных) значений, так же, раскладывая дискретным преобразованием Фурье набор из N точек, мы обязаны получить в аккурат N значений коэффициентов при синусах-косинусах (с небольшими поправками, о которых уже говорилось).
 
А вот вопрос про  
 
Цитата:
могу ли я (как?), имея набор n дискретных данных Фi в диапазоне T1 .. Tn
при помощи DFT посчитать _наперёд заданное количество_ m Fj точек в диапазоне W1 .. Wm

- уже смахивает на задачу интерполяции. Безусловно, можете. 1) Раскладываете свои данные в ряд Фурье - ну, фактически, в тригонометрический полином; 2) считаете значения этого полинома во всех интересующих точках.  
 
Ограничение ряда Фурье одними синусами или косинусами (как у Вас) (в зависимости от "четности" функции) - обычная практика. Коэффициентов тогда получится, понятно, N/2. Главное, чтобы поведение функции соответствовало этой модели. Т.е. если хотим записать функцию косинусами - она должна начинаться с нуля и заканчиваться в ноль.

Всего записей: 2485 | Зарегистр. 21-06-2005 | Отправлено: 14:42 11-04-2011
akaGM

Platinum Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
XPEHOMETP
 
хорошо, спасибо...
 
как в старом добром анекдоте: "пошёл слушать свои валенки" == считать как считал, непрерывным интегралом...
 
последний вопрос:
что мне лучше всего подойдёт из fftpack, RFFTF, DFFTF?

Всего записей: 24766 | Зарегистр. 06-12-2002 | Отправлено: 15:03 11-04-2011
XPEHOMETP

Silver Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
akaGM

Цитата:
что мне лучше всего подойдёт из fftpack, RFFTF, DFFTF?

Честно говоря, я не врубился, что именно они выдают в качестве конечного значения. То есть, в комментариях говорится, что они на-гора выдают некий полупродукт, вместо самих коэффициентов синусов-косинусов. Сравнивал с коэффициентами из книги Уиттекера и Робинсона (о чем говорил несколько раньше), ну ни фига не понял, что там при разложении на выходе. Чтобы понять, наверняка надо читать статьи авторов алгоритмов, вроде Paul Swarztrauber. Ссылки есть, но у нашего универа нет электронной подписки на такие старые года. Так что остался в неведении. Конечно, уже упомянутая ezfftf из этого полупродукта как-то делает реальные коэффициенты, но я поленился разбирать, как и что там. Просто взял эту прогу и использую.
 
Хотя, если верить тестам, которые провел John Burkardt, при разложении в ряд Фурье и обратной свертке эти программы дурят куда меньше, чем ezfftf. Та при нечетном числе точек может реально заглючить.
 
ЗЫ: Кажется, постиг причину взаимного неполного понимания. Ну, это еще Vskazka писал:

Цитата:
 Преобразование Фурье (просто, как таковое) - одно.

Чем-то подобным (преобразованием в общем виде, типа, бесконечным рядом, в пределе - интегралом) Вы и пользовались. Заменяя интеграл некоей суммой по произвольному количеству точек, по своему усмотрению. Имели право. Но я в эту конкретную область рядов Фурье просто не влезал. Мне для практического использования дискретного преобразования Фурье с головой хватает.

Всего записей: 2485 | Зарегистр. 21-06-2005 | Отправлено: 15:31 11-04-2011 | Исправлено: XPEHOMETP, 16:21 11-04-2011
akaGM

Platinum Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
пощупал FFTpack (CFFTF/CFFTB), библиотека хорошая, осталось только расставить везде real*8 :)
 
одно но:
мало того, что я теперь ограничен размером расчитываемых данных (= входным), так ведь ещё половина счёта -- это зеркальное отражение нижних частот, т.е. теперь вообще
N_расчётное = N/2 входное :(

Всего записей: 24766 | Зарегистр. 06-12-2002 | Отправлено: 23:42 11-04-2011
Andrew10

Advanced Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
akaGM
Тебе принципиально иметь исходный код? Если нет, то почему бы ни воспользоваться модулями, входящими в MKL или IMSL?

Всего записей: 780 | Зарегистр. 26-02-2005 | Отправлено: 01:02 12-04-2011
akaGM

Platinum Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
Andrew10

Цитата:
MKL или IMSL?

они ж коммерческие, а мои вещи (пусть звучит нескромно) по миру ходят :)
да и не люблю я их... потому как не в исходниках они :)
 
XPEHOMETP

Цитата:
я уже сижу и на real*8 переправляю

Цитата:
Но - все одинарной точности. Причем указать при компиляции ключиком, что real = real*8, не прокатит: там внутри постоянные прописаны через Е

оказалось, что можно и не править, т.е. прокатило:
ifort  /real-size:64 fftpack.f90

Всего записей: 24766 | Зарегистр. 06-12-2002 | Отправлено: 01:11 12-04-2011 | Исправлено: akaGM, 03:06 12-04-2011
Andrew10

Advanced Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
akaGM

Цитата:
они ж коммерческие, а мои вещи (пусть звучит нескромно) по миру ходят
да и не люблю я их... потому как не в исходниках они  

 
Насколько я помню, в MKL для ДПФ используетcя код библиотеки FFTW-3, которая распространяется в исходниках под лицензией GPL.

Всего записей: 780 | Зарегистр. 26-02-2005 | Отправлено: 12:06 12-04-2011
XPEHOMETP

Silver Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
Andrew10

Цитата:
Насколько я помню, в MKL для ДПФ используетcя код библиотеки FFTW-3, которая распространяется в исходниках под лицензией GPL.

FFTW - они же извращенцы, считают через экспоненты с мнимым показателем... Настоящие пацаны считают через синусы-косинусы, потому что это быстрее! Да, реально: если найдешь реализацию FFT со счетом через экспоненты, то это точно будет под С++ или Дельфи, которые и так тормоз по сравнению с Фортраном...

Всего записей: 2485 | Зарегистр. 21-06-2005 | Отправлено: 13:04 12-04-2011
akaGM

Platinum Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
Andrew10

Цитата:
Насколько я помню, в MKL для ДПФ используетcя код библиотеки FFTW-3

ну так и на фига козе баян?
если "лучшая на западе" в сорсах лежит...
 
http://www.fftw.org/
 
---
по поводу sin/cos ХРЕНОМЕТР прав:
 
у fpu инструции fsin, fcos (и даже fsincos) есть, а fexp -- нет...

Всего записей: 24766 | Зарегистр. 06-12-2002 | Отправлено: 13:16 12-04-2011 | Исправлено: akaGM, 13:35 12-04-2011
terminat0r



Silver Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
akaGM

Цитата:
у fpu инструции fsin, fcos (и даже fsincos) есть, а fexp -- нет...

Вы там не очень розганяйтесь на интеловские fsin, fcos, они точны до 10 знаков только и по сей день. Старый интеловский баг

Всего записей: 2084 | Зарегистр. 31-03-2002 | Отправлено: 20:46 12-04-2011
akaGM

Platinum Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
можно подумать на ителовских камнях есть другая тригонометрия...

Всего записей: 24766 | Зарегистр. 06-12-2002 | Отправлено: 21:32 12-04-2011
terminat0r



Silver Member
Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
akaGM

Цитата:
можно подумать на ителовских камнях есть другая тригонометрия...

нет, но все "умные библиотеки" и компиляторы уже как 20 лет не считают синус с помощью сопроцессора
 
Добавлено:
Вот такой фокус гуляет давно в интернетах

Цитата:
#include <stdio.h>
#include <math.h>
double fsin( double x )
{
register double res;
asm( "fsin" : "=t" (res) : "0" (x) );
return res;
}
 
int
main(int argc, char** argv)
{
printf("sin(43998769152.000000) is %.16le\n", sin(43998769152.000000));
printf("fsin(43998769152.000000) is %.16le\n", fsin(43998769152.000000));
}

 
term@box:~> gcc test.c -o test.x
term@box:~> ./test.x  
sin(43998769152.000000) is -4.0252920638371055e-09
fsin(43998769152.000000) is -4.0819367251002277e-09
 
term@box:~> icc test.c -o test.x
term@box:~> ./test.x
sin(43998769152.000000) is -4.0252920638371055e-09
fsin(43998769152.000000) is -4.0819367251002277e-09

Всего записей: 2084 | Зарегистр. 31-03-2002 | Отправлено: 00:26 13-04-2011 | Исправлено: terminat0r, 00:29 13-04-2011
Открыть новую тему     Написать ответ в эту тему

Страницы: 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 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330

Компьютерный форум Ru.Board » Компьютеры » Прикладное программирование » Закладки » Вопросы программирования на FORTRAN (ФОРТРАН)


Реклама на форуме Ru.Board.

Powered by Ikonboard "v2.1.7b" © 2000 Ikonboard.com
Modified by Ru.B0ard
© Ru.B0ard 2000-2024

BitCoin: 1NGG1chHtUvrtEqjeerQCKDMUi6S6CG4iC

Рейтинг.ru