помогите с ДПФ

[lukashev_sergey]

Какая связь между дискретным преобразованием Фурье и непрерывным? Можно ли получить в matlab тот же результат, что и на бумаге? У меня что-то не выходит. Как нужно нормировать по амплитуде то, что выдает fft? В help-е к matlab-y они предлагают делить на количество ненулевых отсчетов:

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
       see here  L = 1000;           % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); 
y = x + 2*randn(size(t));     % Sinusoids plus noise
plot(Fs*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
        see here    Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

Для синусов это действительно дает верный результат, но , например, для ступеньки это уже не так. Кто-нибудь может привести правильный нормировочный множитель и показать почему он должен быть такой.
Заранее спасибо.

[profrotter]

Преобразование Фурье непрерывного сигнала $s(t)$ определяется выражением:

$$S(\omega)=\int\limits_{-\infty}^{+\infty}s(t)e^{-i\omega t}dt.$$

Функция $S(\omega)$ называется спектральной плотностью сигнала. Преобразование Фурье может быть получено также и для дискртеного сигнала $s_d(t)=\sum\limits_{n=0}^{N-1}s(nT)\delta(t-nT)$, где $s(nT)$ - отсчёты непрерывного сигнала, взятые с периодом $T\leqslant \frac {2\pi} {\omega_m}, \omega_m$ - максимальная частота в спектре сигнала, $N$ - количество ненулевых отсчётов:

$$S_d(\omega)=\sum\limits_{n=0}^{N-1}s(nT)e^{-i\omega nT}.$$

Отсчёты спектра дискретного сигнала следует брать с шагом $\Omega\leqslant \frac {2\pi} {T_c}$, где $T_c=(N-1)T$ - длительность сигнала. Отсчёт спектра дискретного сигнала тогда даётся формулой:

$$S_d(k\Omega)=\sum\limits_{n=0}^{N-1}s(nT)e^{-ikn\Omega T}.$$

Рассмотрим отсчёты спектра дискретного сигнала, взятые с максимальным шагом $\Omega=\frac {2\pi} {T_c}$.При большом количестве отсчётов $N-1\approx N;T_c=NT$ и $\Omega T=\frac {2\pi} N$, тогда

$$S_d(k\Omega)=\sum\limits_{n=0}^{N-1}s(nT)e^{-\frac {2\pi i} N kn}.$$

Полученное выражение называется дискретным преобразованием Фурье (ДПФ) и даёт отсчёты спектра дискретного сигнала, взятые с максимальным шагом.
Спектры дискретного и непрерывного сигналов между собой связаны: спектральная плотность дискретного сигнала представляет собой непрерывную функцию, образованную периодическим повторением спектра непрерывного сигнала:

$$S_d(\omega)=\frac 1 T\sum\limits_{k=-\infty}^{+\infty}S(\omega -k\frac {2\pi} T).$$

с периодом по частоте $\frac {2\pi} T$. Если дискретизация сигнала выполнена корректно (то есть $T\leqslant \frac {2\pi} {\omega_m})$, то можно считать, что на интервале частот $-\frac \pi T\leqslant\omega\leqslant\frac \pi T$ спектр дискретного сигнала пропорционален спектру непрерывного сигнала:

$$S(\omega)\approx T S_d(\omega),$$

соответственно отсчёты спектра дискретного сигнала, взятые с максимальным шагом (ДПФ) пропорциональны отсчётам спектра непрерывного сигнала:

$$S(k\Omega)\approx T S_d(k\Omega).$$

Подробнее:
1. А.Н. Денисенко Сигналы. (примерно 2003г.)
2. Рабинер, Гоулд Теория и применение цифровой обработки сигналов. (примерно 1980г.)
3. Карташов В.Г. Основы теории дискретных сигналов и цифровых фильтров. (примерно 1980г.)

[lukashev_sergey]

Большое списибо, profrotter.

Цитата:

$$S_d(\omega)=\frac 1 T\sum\limits_{k=-\infty}^{+\infty}S(\omega -k\frac {2\pi} T).$$

Это равенство многое объясняет. Стыдно признаться, но оно мне встречалось, и не раз. Я уже методом тыка установил, что для ступеньки множитель $T$ действительно дает верный результат. А как объяснить нормировочный коэффициент когда сигнал не интегрируемый, синус например?

[chainreaction]

Есле владеете английским, то неплохая серия вводных лекций по ДПФ в Matlab есть тут
http://blinkdagger.com/matlab/matlab-introductory-fft-tutorial/

Объясняется на примерах нормировка, отличие одностороннего от двухстороннего спектра и как это реализовать. Даётся пример функции, которая на выходе даёт "интуитивный" спектр.

Вот ссылки на все лекции из этой серии:
http://blinkdagger.com/category/control-systems/

[profrotter]

lukashev_sergey в сообщении #428251 писал(а):

Это равенство многое объясняет. Стыдно признаться, но оно мне встречалось, и не раз. Я уже методом тыка установил, что для ступеньки множитель действительно дает верный результат. А как объяснить нормировочный коэффициент когда сигнал не интегрируемый, синус например?

Ступенька тоже не интегрируемая, если что. Надеюсь вы уверены, что вас интересуют именно спектральные плотности этих сигналов.
Если вы рассматриваете "ступеньку" $\sigma(t)$ или синус $s(t)=\sigma(t)sin(\omega_0 t)$, то для них сначала следует определить преобразование Лапласа, затем можно перейти к преобразованию Фурье, с учётом того, что подынтегральная функция обратного преобразования Лапласа будет иметь полюс на контуре интегрирования. Это приведёт к появлению дельта-функции в выражении для спектральной плотности. (см. И.С. Гоноровский Радиотехнические цепи и сигналы, 1986, приложение I; Свешников, Тихонов Теория функций комплексного переменного). Спектральная плотность ступенчатой функции получится в виде: $S_{\sigma}(\omega)=\pi \delta(\omega)+\frac 1 {i\omega}$. Спектральную плотность синуса можно получить из предыдущего выражения с учётом теоремы смещения спектра или опять же переходя от Лапласова изображения: $S(\omega)=-i\frac {\pi} 2 (\delta(\omega-\omega_0)-\delta(\omega+\omega_0))+\frac {\omega_0} {\omega_0^2 - \omega^2}$.
В случае дискретной ступенчатой функции или одностороннего синуса следует отыскать Z - преобразование, после чего перейти к преобразованию Фурье, учитывая опять же, что такой преход означает, что интегрирование в обратном Z - преобразовании будет происходить по границе круга, вне которого оно сходится.
Результат для ступенчатой функции $S_d(\omega)=\sum\limits_{k=-\infty}^{+\infty}\delta(\omega - \frac {2\pi k} T)+\frac 1 {1-e^{-i\omega T}}$. Для синуса, если есть необходиомость, получите сами. Найдёте выражение для Z - преобразования дискретного синуса $s(nT)=\sigma(nT)sin(\omega_0 nT)$, в нём сделаете замену $z=e^{i\omega T}$, дельта функции будут появляться на частотах $\frac {2\pi k} T \pm \omega_0, k=0,1,2...$
Собственно при всём при этом связь между спектрами аналогового и дискретного сигналов остаётся в силе, соответственно и не изменится нормировочный множитель. Дискретизация спектра для рассматриваемых случаев не может быть корректной , ибо они касаются сигналов бесконечной длительности, соответственно не имеет смысла и ДПФ, ибо оно рассматриватеся для конечного числа отсчётов. Интересно, как вы заставите матлаб считать ДПФ при бесконечном числе отсчётов?
Попробуйте к нормировочному множителю подойти с другой стороны. Представьте, что мы хотим посчитать интеграл Фурье численно методом прямоугольников, тогда мы запишем:

$$S(\omega)\approx T\sum\limits_{n=0}^{N-1}s(nT)e^{-i\omega nT}.$$

Сравните это с выражением для $S_d(\omega)$, которое я приводил в первом сообщении.

[lukashev_sergey]

Напрасно я употребил слово "ступенька" . Вы, видимо, подумали что речь идет о функции Хэвисайда, а я имел в виду функцию вида $1, x\in[0,L]$ и $0$ -- в противном случае. Для этой функции если нормировать fft умножением на $T$, то мы увидим ее преобразование Фурье $ L e^{-\frac{i \omega L }{2}} \frac{ \sin{\frac {\omega L }{2}}} {\frac{\omega L }{2} } $ . С другой стороны, если взять синус и нормировать его fft делением на число ненулевых отсчетов, то в результате мы увидим его разложение в ряд Фурье (домножение на $ T$ не поможет), что везде интерпретируется как правильный результат. Я хочу знать откуда следует такая нормировка. Правильно ли я понимаю что разница между этими двумя сигналами в том, что синус неинтегрируемый? Его спектральная плотность с точки зрения преобразования Фурье это две дельта функции.

[profrotter]

lukashev_sergey в сообщении #428583 писал(а):

... функцию вида $1, x\in[0,L]$ и $0$ -- в противном случае. Для этой функции если нормировать fft умножением на $T$, то мы увидим ее преобразование Фурье $ L e^{-\frac{i \omega L }{2}} \frac{ \sin{\frac {\omega L }{2}}} {\frac{\omega L }{2} } $ .

Тут можно всё увидеть наглядно. Спектральная плотность для дискретного прямоугольного импульса из $N$ отсчётов даётся выражением

$$S_d(\omega)=\frac {sin(\frac {N\omega T} 2)} {sin(\frac {\omega T} 2)} e^{-i\frac {\omega N T} 2}.$$

В нуле амплитудный спектр принимает значение $N$. Для аналогового прямоугольного импульса длительностью $L$ спектральная плотность определена выражением, которое вы привели. В нуле получается значение $L$. Теперь вспоминаем, что $L\approx NT$. Можете построить графики $S(\omega)$ и $TS_d(\omega)$ и посмотреть на них.

Тут красным спектр непрерывного сигнала, синим - дискретного, зелёным - результат ДПФ, 20 отсчётов в имульсе. Красный и синий в нуле не совпадают, потому что на самом деле $L=(N-1)T$ и при 20-ти отсчётах это имеет значение.

lukashev_sergey в сообщении #428583 писал(а):

С другой стороны, если взять синус и нормировать его fft делением на число ненулевых отсчетов, то в результате мы увидим его разложение в ряд Фурье (домножение на $ T$ не поможет), что везде интерпретируется как правильный результат. Я хочу знать откуда следует такая нормировка. Правильно ли я понимаю что разница между этими двумя сигналами в том, что синус неинтегрируемый? Его спектральная плотность с точки зрения преобразования Фурье это две дельта функции.

Периодический сигнал с периодом $T_{\Omega}=\frac {2\pi} {\Omega}$ раскладывается в комплексный ряд Фурье:

$$s_p(t)=\sum\limits_{n=-\infty}^{+\infty}C_ne^{i n \Omega t}.$$

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

$$C_k=\frac 1 {T_{\Omega}} S(k\Omega).$$

То есть коэффициенты ряда Фурье периодического сигнала можно получить дискретизацией спектра соответствующего непериодического сигнала с шагом, равным частоте повторения. А теперь так:

$$C_k=\frac 1 {T_{\Omega}} S(k\Omega)\approx \frac T {T_{\Omega}} S_d(k\Omega).$$

При ДПФ $\Omega=\frac {2\pi} L;T_{\Omega}=\frac {2\pi} {\Omega}=L\approx NT$, откуда

$$C_k\approx \frac T {NT} S_d(k\Omega)=\frac 1 N S_d(k\Omega).$$

Под $S_d(k\Omega)$ тут следует понимать ДПФ отсчётов периодического сигнала на периоде.

[lukashev_sergey]

Спасибо за объяснение, profrotter.
В начале укажу замеченные неточности:

profrotter в сообщении #428712 писал(а):

Спектральная плотность для дискретного прямоугольного импульса из $N$ отсчётов даётся выражением

$$S_d(\omega)=\frac {sin(\frac {N\omega T} 2)} {sin(\frac {\omega T} 2)} e^{-i\frac {\omega N T} 2}.$$

Помоему, если первый отчет соответствует $t=0$, то это выражение должно выглядеть так

$$S_d(\omega)=\frac {sin(\frac {N\omega T} 2)} {sin(\frac {\omega T} 2)} e^{-i\frac {\omega (N-1) T} 2}$$

profrotter в сообщении #428712 писал(а):

Тут красным спектр непрерывного сигнала, синим - дискретного, зелёным - результат ДПФ, 20 отсчётов в имульсе.

Что-то зеленого не видно .

Вы полностью ответили на мой вопрос.

profrotter в сообщении #428712 писал(а):

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

$$C_k=\frac 1 {T_{\Omega}} S(k\Omega).$$

Выражение$S(\omega)=\int\limits_{-\infty}^{+\infty}s(t)e^{-i\omega t}dt$ при $\omega=\frac{2\pi}{T_c}k$ и соответствующих пределах с точностью до множителя, нормирующего экспоненту, действительно соответствует формуле для коэффициентов ряда Фурье. Теперь, кстати, понятно что имеется в виду, когда говорят что ДПФ соответствует периодическому продолжению исходного сигнала: последнее утверждение используется в обратном направлении.
Помоему, в моей голове сейчас сложилось верное представление о том, что такое ДПФ и как его интерпретировать.

[profrotter]

lukashev_sergey в сообщении #428811 писал(а):

Помоему, если первый отчет соответствует , то это выражение должно выглядеть так

Да.

lukashev_sergey в сообщении #428811 писал(а):

Что-то зеленого не видно.

Зелёный - это точки. Ищите их на оси абсцисс, а на частотах 0, -120, 120 изображены выколотые "точки с подставкой". Впрочем такой же график вы можете построить в матлаб.
Думаю разобрались.

[lukashev_sergey]

А! Понял! Число шагов подобрано так, что отсчеты с нулями совпали.