Точно измерить частоту. Опции

[blackfin]

Просто здаётся мне: что БПФ, что автокореляция - это много процессорного времени.

Хмм.. много или не много?.. Complex 16-bit FFT:

Код

Performance for BF527:
Cycle Count:
                  10008091 cycles for FFT size of  8192. (data in SDRAM + exp in SDRAM)
                  21268141 cycles for FFT size of 16384. (data in SDRAM + exp in SDRAM)
                  44727982 cycles for FFT size of 32768. (data in SDRAM + exp in SDRAM)
                 102921037 cycles for FFT size of 65536. (data in SDRAM + exp in SDRAM)

Performance for BF547:
Cycle Count:
                   5208301 cycles for FFT size of  8192. (data in SDRAM + exp in SDRAM)
                  10969341 cycles for FFT size of 16384. (data in SDRAM + exp in SDRAM)
                  22768076 cycles for FFT size of 32768. (data in SDRAM + exp in SDRAM)
                  53398332 cycles for FFT size of 65536. (data in SDRAM + exp in SDRAM)

[bahurin]

Такая обработка ничем не лучше БПФ, хотя по определению должна давать точность по максимуму правдоподобия (CRLB).

Интересно, почему тогда системы слежения выделяют фазу и производят подстройку частоты именно на основе фазового измерения, а не на основе БПФ?
На самом деле такая оценка дает точность соответсвующую методу максимального правдоподобия. Именно поэтому она лучше БПФ. При бпф после умножения на комплексную экспоненту все значения складываются (этакий фильтр) в результате вы получите некую среднюю фазу и не сможете понять как она меняется на интервале обработки. Здесь же у вас есть информация о фазе. Зная что у нас одна синусоида это линейная функция. Это можно использовать для определения df. Кстати в этом случае раскрыть периодичность арктангенса вполне не сложно читай например здесь

ЗЫ кстати данная реализация куда быстее бпф, так как арктангенс можно вообще не вычислять а задать таблично.

Сообщение отредактировал bahurin - Mar 3 2010, 08:41

[fontp]

Интересно, почему тогда системы слежения выделяют фазу и производят подстройку частоты именно на основе фазового измерения, а не на основе БПФ?
На самом деле такая оценка дает точность соответсвующую методу максимального правдоподобия. Именно поэтому она лучше БПФ. При бпф после умножения на комплексную экспоненту все значения складываются (этакий фильтр) в результате вы получите некую среднюю фазу и не сможете понять как она меняется на интервале обработки. Здесь же у вас есть информация о фазе. Зная что у нас одна синусоида это линейная функция. Это можно использовать для определения df. Кстати в этом случае раскрыть периодичность арктангенса вполне не сложно читай например здесь

ЗЫ кстати данная реализация куда быстее бпф, так как арктангенс можно вообще не вычислять а задать таблично.

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

Возможно "такая оценка" с измерением фазы и будет лучше, но только при ещё более высоких SNR.
При SNR<10 дб "такая оценка" наоборот работать будет плохо, поскольку сама возможность восстанавливать линейную фазу по непрерывности по значению в пределах 2*pi исчезает - всё забито шумом. Это я говорю по собственному опыту. Я все эти способы пробовал.

По производительности. Это арктангенс двух переменных с восстановлением непрерывности фазы
Ph[m] = arctg2(SIm,SQm);
if (ABS(Ph[m-1]-Ph[m]) > M_PI) {// phase jump is detected
if (ABS(Ph[m-1]-Ph[m]-2*M_PI) <= M_PI) {
Ph[m] += 2*M_PI;
}else {
Ph[m] -= 2*M_PI;
}
}

1. 2 сравнения чтобы выяснить квадрант
2. Деление (x/y) или (y/x) после сравнения кто из них меньше, чтобы не делить на примерно 0 и не получать примерно бесконечность
3. Чтобы считать арктангенс уже от одной переменной в лоб по таблице с приличной точностью таблица должна быть уж очень большой. Поэтому обычно считают по быстрому алгоритму CORDIC по небольшой "фрактальной" таблице для одного квадранта.
4. Восстановление непрерывности фазы - еще несколько сравнений и вычитаний как написано выше

Всё про всё займет если не аппаратно тактов 50-100. Итого имеем пусть 50*N - N длина блока. Вычисление FFT требует примерно 4*N*log(N), DFT на современной архитектуре с двумя МАС-аккумуляторами - просто k*N, k - число вычисляемых гармоник в диапазоне. Плюс корреляция (Ph[i]*i), необходимая для подгонки линейной модели (регрессии, сумма(i*i) пусть расчитано предварительно, частота по линейной регрессии = сумма(Ph[i]*i)/сумма(i*i)). Арктангенсы выгодней считать получается только ассимптотически при очень больших блоках N и то только в некоторых случаях, когда 50 гармоник DFT (синхронных детекторов) недостаточно.

ЗЫ. Что касается "систем слежения" (видимо ФАПЧ c обратной связью) то можно совершенно определённо сказать, что всякий ФАПЧ в синхронизм входит медленней, чем это делают измерительные feed-forward системы, причем последние для области низких SNR работают как раз чаще всего на оценивании спектра через интерполяцию.
Если не ставить целью провести измерение частоты максимально быстро, то в принципе, задачу можно решать стандартно именно с помощью ФАПЧ. Вполне приемлемое решение, но оно никогда не будет оптимально по времени измерения. Есть даже такие готовые микросхемы квадратурной демодуляции с фазовым детектором и цифровым генератором. Код управляющий частотой и соответствует частоте входного сигнала. Только это всегда будет немного тормозное решение, уступающее схемам по прямому измерению без обратной связи по точности/времени измерения

[petrov]

Она не лучше спектрального оценивания ("БПФ"), именно потому, что в очень широких пределах SNR (0-40дб) после интерполяции спектра это
1. 2 сравнения чтобы выяснить квадрант
2. Деление (x/y) или (y/x) после сравнения кто из них меньше, бо тот кто их них больше уйдет в разнос
3. Чтобы считать арктангенс уже одной переменной в лоб по таблице с приличной точностью таблица должна быть очень большой. Поэтому обычно считают по быстрому алгоритму CORDIC по небольшой "фрактальной" таблице для одного квадранта.

Без всякого деления и арктангенса с одним сравнением сразу кордиком посчитаем аргумент комплексного числа.

Всё про всё займет если не аппаратно тактов 50-100. Итого имеем 50*N - N длина блока.

Прежде чем "арктангенсы" считать ещё проредим простейшим скользящим средним CIC-ом раз в 10.

[fontp]

Без всякого деления и арктангенса с одним сравнением сразу кордиком посчитаем аргумент комплексного числа.
Прежде чем "арктангенсы" считать ещё проредим простейшим скользящим средним CIC-ом раз в 10.

Кордик в любом случае даст десятки тактов. Скорее 50, чем 10, поскольку там тоже всё на сравнениях. А для DSP сравнения значительно затратнее МАС. Я просто пример привёл с СORDIC-ом на один квадрант на вскидку.
CIC может решать. Раз в десять только, если диапазон измеряемых частот в десять раз меньше частоты Найквиста.
Там даже не CIC, а простое усреднение нужно до вычисления арктангенса, если блочный алгоритм. Усреднять возможно только пока не подавляется синком полезный сигнал. Усреднение ещё и полезно, чтобы не позволить шуму делать сбои при восстановлении фазы по непрерывности, увеличивает SNR. Усреднение принято ))
Но похожую предварительную обработку можно проделать и для DFT, это не аргумент при сравнении.

Хорошо. Такой алгоритм будет оптимален по точности и вычислительно. Но не такой робастный как спектральное оценивание и грохнется раньше при снижении SNR. Из-за необходимости устранять скачки фазы

[petrov]

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

Но именно так в модемах и делают, работает при плохом сигнал/шум, грохнуться ведь в принципе любой метод может, в общем прям так сразу не очевидно что спектральное оценивание робастней.

[fontp]

Но именно так в модемах и делают, работает при плохом сигнал/шум, грохнуться ведь в принципе любой метод может, в общем прям так сразу не очевидно что спектральное оценивание робастней.

В BURST-модемах делают все три метода (спектр, корреляция и линейная фаза). В общем-то все три уязвимы.
Но особенно тот, который строит модель линейной фазы по циклической.
Шумовой импульс по фазе в любой точке с амплитудой порядка ПИ погубит модель.
Вам известен простой алгоритм с интегральными свойствами? Чтобы строил линейную модель игнорируя одиночные выбросы. Так сказать алгоритм циклической регрессии, но чтобы без перебора. Мне нет.
(Хотя и можно, наверное, что то такое приладить как вычисление разности фаз между соседними точками с последующим ранговым усреднением медианой. Но это уже будет много более затратно)

Восстановление фазы по непрерывности я пробовал и мне сравнительно не понравилось. Очень быстро и резко кончается.
В спектре Фурье одиночные импульсы размазываются и никому не мешают. Хотя, конечно, наоборот гармониченская помеха в рабочем диапазоне критична.

[bahurin]

Возможно "такая оценка" с измерением фазы и будет лучше, но только при ещё более высоких SNR.
При SNR<10 дб "такая оценка" наоборот работать будет плохо, поскольку сама возможность восстанавливать линейную фазу по непрерывности по значению в пределах 2*pi исчезает - всё забито шумом. Это я говорю по собственному опыту. Я все эти способы пробовал.

Думаю, что поставленные требования к точности оценки частоты подразумевают, что отношение сигнал-шум велико (гораздо больше 10 дБ). Именно поэтому я и предложил такой вариант.

Чтобы считать арктангенс уже от одной переменной в лоб по таблице с приличной точностью таблица должна быть уж очень большой

Как раз в этом случае высокая точность арктангенса и не нужна, поскольку мы измеряем не мгновенную фазу, а "медленное" изменение этой фазы. Так например если есть таблица из 360-и значений арктангенса, то средняя ошибка измерения одной фазы будет 0.5 градуса. 0.5 градуса на интервале T = 1 сек приведут к ошибке по частоте равной 0.5*pi/180 =0.008 Гц.

[fontp]

Думаю, что поставленные требования к точности оценки частоты подразумевают, что отношение сигнал-шум велико (гораздо больше 10 дБ). Именно поэтому я и предложил такой вариант.

Для высокого отношения сигнал/шум метод может оказаться хорошим

Как раз в этом случае высокая точность арктангенса и не нужна, поскольку мы измеряем не мгновенную фазу, а "медленное" изменение этой фазы. Так например если есть таблица из 360-и значений арктангенса, то средняя ошибка измерения одной фазы будет 0.5 градуса. 0.5 градуса на интервале T = 1 сек приведут к ошибке по частоте равной 0.5*pi/180 =0.008 Гц.

1 сек это не по максимуму правдоподобия.
По максимуму правдоподобия это 1/8 сек (у TigerSHARCa), согласно критерию при 20 дб, а при большем отношении SNR ещё меньше.
У Вас получится примерно в среднем sqrt(2)*0.008*8 (sqrt(2) поскольку разность фаз на интервале по Вашей логике по двум точкам в конце и начале интервала), т.е. за 1/8 секунды по двум арктангенсам Вы так получите почти 0.1 гц ошибки. Мы ведь говорим о предельных оценках. Теория говорит, что
var(F) = 3*(Fs*Fs)/(2*pi*pi*N*(N-1)*(N-1)*(Es/No))
При этом нужная таблица получится в 10 раз больше для точности 0.01 гц если по двум точкам считать.
(По двум точкам на границах всего интервала считать нельзя, поскольку набег фазы должен быть меньше pi,
но для оценок, ладно сойдет. Мы будем на самом деле считать разность фаз между соседними арктангенсами, а потом усредним по всем разностям)

На самом деле строя линейную регрессию по минимуму квадратов по многим измерениям фазы точность повышается примерно как корень из числа измерений и по отношению к систематическим ошибкам в таблице. Здесь есть ещё небольшой резерв в несколько раз. В любом случае таблица нужна в несколько раз больше. И всё это ещё нужно очень хорошо сбалансировать, по равенству (<=) систематической ошибки примерно случайной, вызваной физическим шумом - посредством баланса по сколько отсчетов усреднять и сколько арктангенсов брать для регрессии=(полное число отсчетов)/(размер лаптя усреднения до вычисления фазы).

[leksa]

Всем здравствуйте!
Тема интересная, позвольте выложить матлабовский скрипт, в котором я набросал алгоритм оценки частоты комплексной экспоненты на основе планарной фильтрации.
Ничего сверхестественного, просто один из adhoc способов. Идея в том, чтобы умножая отсчеты комплексной экспоненты на свою задержанную и взятую с комплексным сопряжением копию, получить последовательность фазоров, угол которых пропорционален частоте исходной экспоненты. Затем нужно сложить полученные фазоры как комплексные числа, получив на комплексной плоскости длиный вектор, угол которого опять же пропорционален искомой частоте(эта операция и есть планарная фильтрация). Таким образом уменьшается воздействие шума на оценку. Далее вычисляем угол этого длинного вектора, умножаем на коэффициент Fs/(2*pi) и получаем оценку частоты.
Правда полученная таким образом оценка достаточно сильно подвержена влиянию шума.
Чтобы несколько улучшить оценку, я пошел немного дальше и используя вычисленную оценку частоты, создаю комплексный КИХ фильтр с ИХ в виде 1 периода комплексной экспоненты на полученной оценкой частоте. Получается такой комплексный ППФ со средней частотой совпадающей с нашей оценкой частоты.
Затем пропускаю входной сигнал через этот ППФ и снова нахожу оценку по вышеописанному методу.
Засчет второй итерации после фильтрации, погрешность оценки может уменьшиться на порядок по сравнению с первой оценкой.
Вот такой вот способ.
Правда, вторая итерация имеет смысл только в определенном диапазоне отношения сигнал шум (если ОСШ выше 45 дБ, то полезный эффект от второй итерации уже незначителен)
Планарная фильтрация используется в ряде прямых(feedforward) методов синхронизации, описанных например в
в часто упоминающейся на форуме книге Synchronization Techniques for Digital Receivers.

Собственно, скрипт:

%test_carr_freq_ff_est
script
clear
clc
close all

Fs=6400;
Tau=1.00;%
Nlen=Tau*Fs;
t=[1/Fs:1/Fs:Nlen/Fs].';
SNR=6;
Fc=47.333333333333333333333;%
disp(['SNR = ' num2str(SNR)])
disp(['Fc = ' num2str(Fc)])
ini_phase = (rand(1)*2-1)*pi
y = 1.*exp(j*(2*pi*Fc*t+ini_phase));
Sinp = awgn(y,SNR,'measured');
Sinp_d = Sinp(1:end-1).*conj(Sinp(2:end));
Sinp_d_pf = sum(Sinp_d)
% цикл чисто для наглядности работы алгоритма
Sinp_d_pf_k = zeros(size(Sinp_d));
for k=2:length(Sinp_d)
Sinp_d_pf_k(k)=Sinp_d_pf_k(k-1)+Sinp_d(k);
% figure(101)
% plot([Sinp_d_pf_k(k-1) Sinp_d_pf_k(k)],'y- *'),grid on,hold on
% pause(1/25)
end
figure(11)
plot(Sinp_d_pf_k,'y- *'),grid on,hold on,title('планарная фильтрация');
phase_pf = -angle(Sinp_d_pf);
phase_pf_pi = phase_pf/pi;
Fc_est = phase_pf_pi*Fs/2;
delta_F = Fc_est-Fc;
disp(['Оценка частоты: Fc_est = ' num2str(Fc_est)])
disp(['Погрешность оценки частоты: delta_F = ' num2str(delta_F)])
figure(12)
sz=prod(size(Sinp));
tv=[1/Fs:1/Fs:sz/Fs];
plot(tv,real(Sinp),'b- .'),grid on,hold on,title('IQ сигнала после фильтра');
plot(tv,imag(Sinp),'r- .'),grid on,hold off
figure(13)
[PSDPSD,f]=pwelch(Sinp,[],[],[],Fs,'twosided');
plot(f-Fs/2,fftshift(10*log10(PSDPSD/max(PSDPSD))),'y- '),grid on,title('Спектр сигнала после фильтра'),hold on;
% пофильтруем
h_len = round(Fs/Fc_est)*1;
hbpf = exp(j*2*pi*(Fc_est/Fs)*[1:h_len]);
Sinp_f = filter(hbpf,1,Sinp);

% еще разок
Sinp = Sinp_f;
Sinp_d = Sinp(1:end-1).*conj(Sinp(2:end));
Sinp_d_pf = sum(Sinp_d);
% цикл чисто для наглядности работы алгоритма
Sinp_d_pf_k = zeros(size(Sinp_d));
for k=2:length(Sinp_d)
Sinp_d_pf_k(k)=Sinp_d_pf_k(k-1)+Sinp_d(k);
% figure(101)
% plot([Sinp_d_pf_k(k-1) Sinp_d_pf_k(k)],'y- *'),grid on,hold on
% pause(1/25)
end
figure(21)
plot(Sinp_d_pf_k,'y- *'),grid on,hold on,title('планарная фильтрация');
phase_pf = -angle(Sinp_d_pf);
phase_pf_pi = phase_pf/pi;
Fc_est2 = phase_pf_pi*Fs/2;
delta_F2 = Fc_est2-Fc;
disp(['Оценка частоты №2: Fc_est2 = ' num2str(Fc_est2)])
disp(['Погрешность оценки частоты: delta_F2 = ' num2str(delta_F2)])
figure(22)
sz=prod(size(Sinp));
tv=[1/Fs:1/Fs:sz/Fs];
plot(tv,real(Sinp),'b- .'),grid on,hold on,title('IQ входного сигнала');
plot(tv,imag(Sinp),'r- .'),grid on,hold off
figure(23)
[PSDPSD,f]=pwelch(Sinp,[],[],[],Fs,'twosided');
plot(f-Fs/2,fftshift(10*log10(PSDPSD/max(PSDPSD))),'y- '),grid on,title('Спектр вх сигнала'),hold on;
% % характеристики фильтра
% figure(31)
% plot(real(hbpf),'b- .'),grid on,hold on
% plot(imag(hbpf),'r- .'),grid on
% figure(32)
% freqz(hbpf,1)

Сообщение отредактировал leksa - Mar 3 2010, 17:28

[TigerSHARC]

для fontp:
хочу всё таки запустить метод маклеода.

Сигнал - это синусоида частотой 53.33Гц.
Дискретизируем его с частоой 960 Гц.
Берём 128 тоочек. (соответственно частотное разрешене DFT = 7.5Гц)
Делаем DFT - получаем размазаный спектр.

Теперь узнаю максимум - 8 бин DFT (т.е 60Гц)

В качестве вектора из трёх элементов отправляю на Маклеода массив из трёх значений соответсвтующих 7-му, 8-му и 9-му бинам DFT. (как положено - максимум и две точки окрестности)
Маклеод выдаёт -0.02
Стало быть смещение относительно максимума -0.02 бина. Это выходит 60 - 0.02*7,5 = 59,85 Гц

Объясните где я неправ. Может нужно разрешение увеличить?

[fontp]

для fontp:
хочу всё таки запустить метод маклеода.

Сигнал - это синусоида частотой 53.33Гц.
Дискретизируем его с частоой 960 Гц.
Берём 128 тоочек. (соответственно частотное разрешене DFT = 7.5Гц)
Делаем DFT - получаем размазаный спектр.

Теперь узнаю максимум - 8 бин DFT (т.е 60Гц)

В качестве вектора из трёх элементов отправляю на Маклеода массив из трёх значений соответсвтующих 7-му, 8-му и 9-му бинам DFT. (как положено - максимум и две точки окрестности)
Маклеод выдаёт -0.02
Стало быть смещение относительно максимума -0.02 бина. Это выходит 60 - 0.02*7,5 = 59,85 Гц

Объясните где я неправ. Может нужно разрешение увеличить?

Во всех спектральных оценивателях (Маклеоде, Квине) речь идёт о комплексной экспоненте. Непосредственно с действительной синусоидой не будет точно работать, поскольку синусоида - это две комплексных экспоненты. Поэтому для того, чтобы добраться до Маклеода (и всех других) нужно ещё сделать из реальной синусоиды комплексную, ну хотя бы тем же квадратурным детектором. Или подавить влияние второго пика широкими спектральными,например, гауссовыми окнами как описано по ссылке Стенфордского университета. Потом делать комплексное преобразование по типу DFT в пределах [-7.5 7.5]. В любом случае это комплексное DFT, а не действительное.

Средний бин будет соответствовать частотному нулю.
первый бин на стороне положительных частот соответствует нулю, восьмой 52.5

Кстати все три рассмотреных метода с действительной синусоидой работать не будут, не только спектральный. Здесь в теме рассмотрены в разных сообщениях всего 3 работающих почти оптимально метода - интерполяция спектра Фурье вблизи максимума, вычисление набега фазы через arg(сумма по T(S()*Sсопряженное())/T (автокорреляционный) и просто в лоб линейное моделирование набега фазы после квадратичного детектирования. Видите, эти формулы даже не напишешь в форуме, читайте лучше статьи....

Когда-то загружал статью, с вариантом перехода от реальной синусоиды к комплексной
http://electronix.ru/forum/index.php?act=A...st&id=37367

[TigerSHARC]

Во всех спектральных оценивателях (Маклеоде, Квине) речь идёт о комплексной экспоненте. Непосредственно с действительной синусоидой не будет точно работать, поскольку синусоида - это две комплексных экспоненты. Поэтому для того, чтобы добраться до Маклеода (и всех других) нужно ещё сделать из реальной синусоиды комплексную, ну хотя бы тем же квадратурным детектором. Или подавить влияние второго пика широкими спектральными,например, гауссовыми окнами как описано по ссылке Стенфордского университета. Потом делать комплексное преобразование по типу DFT в пределах [-7.5 7.5]. В любом случае это комплексное DFT, а не действительное.

Средний бин будет соответствовать частотному нулю.
первый бин на стороне положительных частот соответствует нулю, восьмой 52.5

Кстати все три рассмотреных метода с действительной синусоидой работать не будут, не только спектральный. Здесь в теме рассмотрены в разных сообщениях всего 3 работающих почти оптимально метода - интерполяция спектра Фурье вблизи максимума, вычисление набега фазы через arg(сумма по T(S()*Sсопряженное())/T (автокорреляционный) и просто в лоб линейное моделирование набега фазы после квадратичного детектирования. Видите, эти формулы даже не напишешь в форуме, читайте лучше статьи....

Когда-то загружал статью, с вариантом перехода от реальной синусоиды к комплексной
http://electronix.ru/forum/index.php?act=A...st&id=37367

Большое спасибо!

[fontp]

вот ещё нашелся тот матлабовский тест по ссылке, "которого больше нет уже в Интернете"
Начните с него с комплексных экспонент, он должен сразу работать
Там три работающих оценщика параметров комплексной экспоненты Macleod, Quinn1, Quinn2

Если пиков в спектре несколько, но они далеко друг от друга - нужно применять окна пусть даже с широким главным максимумом (разрешение нас не волнует), но максимально подавленными боковиками в спектральной области. Это может быть окно Кайзера, экспоненциальное или даже гауссовское. Видимо, длину выборки при применении окон нужно увеличивать по отношению к пороговой. Это работает даже без квадратурного детектирования.

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

Прикрепленные файлы

 tst3.zip ( 4.48 килобайт ) Кол-во скачиваний: 56

 

[blackfin]

Сравните с оцениванием спектра и обратной квадратичной интерполяцией. Если интервал частот узкий, то сразу можно начинать с того, что называют ML-extension

1. Умножить на центральную комплексную экспоненту (снести сигнал в 0) и отфильтровать ФНЧ. Соответственно прорежение.
2. Взять 5-7 сумм типа ДПФ в интервале [-7.5, 7.5]. Оценка энергии
3. Найти максимум энергии и по трём точкам построить параболу. Аргумент максимума - частота

Мне вот, непонятно, зачем "умножить на центральную комплексную экспоненту"? При таком сносе отклонение вправо/влево измеряемого вещественного сигнала от частоты 50Гц станут неотличимы друг от друга.

Мне кажется, что проще "умножить на комплексную экспоненту" с частотой 42.5Гц и т.о., снести в 0 левый участок спектра 42.5Гц после чего уже сделать FFT на интервале [0, 15].

Сообщение отредактировал blackfin - Mar 4 2010, 07:56