Цитата(tocha @ May 3 2011, 12:45)
А теперь хочу сделать это (IRFFT), то есть irfft симметричного комплексного сигнала у которого будет реальный выход. Максимально эффективно с точки зрения вычислительных затрат.
Если можете ткнуть в готовый код, буду благодарен (в пятницу навскидку не нашёл - может плохо искал...).
То есть, если я правильно понял, вы такое (IRFFT) делали: обратное преобразование симметричного сигнала длиной 2*N с помощью N-point ICFFT плюс немного вычислений и этим самым экономили ресурсы? Если скажете да, буду тоже благодарен.
IRFFT как бы прилагается к RFTP в качестве обратной функции, а потому желательно, чтобы они были "от одного производителя". Лучше поищите вы в том же месте, где нашли прямую RFFT.
Тут есть один нюанс: если действительный массив, длиной 2*N, обработать RFTP, а к результату этой обработки сразу применить обратное IRFTP, то может получиться два варианта, в зависимости от реализации алгоритма. В первом варианте вы получите назад свой исходный массив с вычислительной точностью. А во втором варианте вы получите его увеличенным по амплитуде в N раз. Вот из-за этого множителя и возникает "несовместимость" прямого и обратного RFFTP у разных реализаторов. Причем это касается не только RFFT, но и простого CFFT - увеличение в N раз у обоих.
А сама проблема сводится к тому, за чей счет подавление этого множителя N относить, если уж решено сделать так, чтобы после обратного преобразования возвращался точно исходный результат. И тут кто во что горазд. Например, в MatLab'е прямое FFT делают, не вводя поправки на множитель, а внутри процедуры обратного IFFT делают деление всех элементов массива на N в начале или конце процедуры (эффект от места делёжки не зависит). Такое решение кажется несправедливым по отношении к IFFT и сильно похоже на подгон результата. Можно было бы разделить коэффициент N и по "справедливости", если делить после каждого преобразования, (как прямого, так и обратного) на квадратный корень из N. Но это противно, т.к. надо вычислять корень.
Я же задалась целью выяснить, как должно быть на самом деле
. Для чего создала на всю длину массива синусоиду на полный период с известной мне амплитудой 1000 единиц, и посмотрела, какие коэффициенты дает прямая FFT. очевидно, что из такой синусоиды по всем канонам должен получиться единственный мнимый коэффициент, равной 1000, во втором по порядку элементе массива. А оно не получилось... Тогда пришлось модифицировать процедуры так: после FFT умножаю на множитель 2/N, а после обратной IFFT делю на 2. Такой прием дает как правильную амплитуду гармоник, так и поглощает множитель при обратном преобразовании, т.к. (2/N)/2 = 1/N.
Очень возможно, что в этом деле я что-то недопонимаю и изобретаю велосипед, но мне удивительно, что люди почему-то над этими вещами не задумываются, а вполне механически применяют процедуры из книг или библиотек функций. Может быть я тут где-то и пронеслась, но вам я советую ОБЯЗАТЕЛЬНО проверить прямое и обратное FFT на тестовом примере из синусоиды (или косинусоиды) с известной вам амплитудой, чтобы проверить FFT-результат на предмет множителя.
Подробнее обсуждать алгоритмы мне с вами сложно, т.к. вы не пишите ни процессор, на котором собираетесь гонять эти процедуры, ни язык, на котором собираетесь их программировать. А мои наработки для Pentium могут оказаться для ваш излишней информацией, т.к. там есть своя специфика. Например, получать значения "опорных" синусов и косинусов на Пентиуме быстрее через тригонометрическую функцию sincos, входящую в набор инструкций сопроцессора (!), которая за один проход дает оба результата - синус и косинус. Это оказалось быстрее, чем получать синусы и косинусы половинных углов по тригонометрической формуле для половинного угла. А, скажем, для DSP-процессора, у которого тригонометрических функций в наборе инструкций нет, это было бы уже плохим решением. А вот в тех случаях, когда N фиксировано, то набор синусов для N/4 (от 0 до 90 градусов) лучше вообще зашить в ПЗУ и не заниматься их вычислением. Иными словами, здесь есть очень много нюансов, которые зависят от сопутствующих обстоятельств.
ifft с симметричным входом, оптимальный алгоритм
Apr 29 2011, 18:15
