Построение спектра Фурье в matlab

marandr

Уже некоторое время пытаюсь разобраться с построением спектра в matlab. Код вроде не сложный, но... для проверки даю на вход гауссиан ( exp(-(t^2)/2) спектр которого должен совпадать с самим сигналом. На выходе получаю сигнал той же формы, но Уже приблизительно в 50 раз (вектор на выходе той же длины, что и исходный). Не подскажете, в чем может быть косяк?
Кусок кода, вычисляющий спектр (u - исходный вектор, длина порядка 1500):
Puu = fft(u);
U=fftshift(Puu);
Puu=U.*conj(U);
Вывожу sqrt(Puu).

stm7543347

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

Nastia-2008

А ты умножал на шаг сетки?
Попробуй так.
Просто FFT в матлабе этого как раз не делает.
Вот что я имею в виду.
delta_f = 1/(N*delta); %шаг в частотной области
%FFT
G = fftshift(fft(fftshift(g * delta;
%IFFT
g = ifftshift(ifft(ifftshift(G * length(G) * delta_f;

stm2383383

В матлабе шикарный хелп с примерами. Пример скрипта, считающий спектр, а также выводящий его в дБ с подписанными осями, имеется там.

marandr

Если имеется ввиду
Fs = 1000; % Sampling frequency
T = 1/Fs; % Sample time
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:50y(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
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
То в слегка доработанном виде
Fs = 1000; % Sampling frequency
T = 1/Fs; % Sample time
[u,txt,raw] = xlsread('5.csv');
L = length(u); % Length of signal
t = (0:L-1)*T; % Time vector
figure(1)
subplot(211)
plot(Fs*t,u)
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(u,NFFT)/L;
y=fftshift(Y);
f = Fs*linspace(0,1,NFFT);
% Plot single-sided amplitude spectrum.
%figure(2);
subplot(212);
plot(f,2*abs(y(1:NFFT
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
он выдает тот же узкий гаусс. Естественно, это я уже пробовал.

marandr

А какая тут норма должна использоваться? Понятно, деление на максимум амплитуды не изменит картинки. Тогда что? Гугл не помог

ulek

А кто сказал, что в частотной координате гауссиана должна быть такой же ширины в пикселях, что и в пространственных координатах?

marandr

Эмм, fft на выходе дает вектор той же длины, что и входной. Поэтому length(G)=N. Понятно, что прямое и обратное преобразование вместе дадут исходный сигнал) Просто домножение на delta при прямом преобразовании не меняет ширины импульса. Или я тебя не понял?

Nastia-2008

короче я тут накидал. Поиграйся с последней частью. Там как раз твой вариант.
L = 5; % numerical window size
N = 32; % number of samples
delta = L / N; % spatial step
delta_f = 1/(N*delta);
%numerical window in spatial and frequency domain
x = (-N/2 : N/2-1) * delta;
f = (-N/2 : N/2-1) * delta_f;
%Gauss function
a = 1;
g = exp(-pi*a*x.^2);
%Gauss FFT
G = fftshift(fft(fftshift(g * delta;
%IFFT
g2 = ifftshift(ifft(ifftshift(G * N * delta_f;
%plot
figure(1)
subplot(3,1,1)
plot(f, abs(G '--r','linewidth',4);
subplot(3,1,2)
plot(x, g, '-g', 'linewidth',3 ); hold on;
%--------------------Your code------------------
Puu = fft(fftshift(g*delta; %you should multiply by delta! try no to :)
U=fftshift(Puu);
Puu=U.*conj(U);
subplot(3,1,3)
plot(f, sqrt(Puu '--b', 'linewidth', 1 );

Nastia-2008

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

marandr

По хорошему a = 1/(2*pi); тогда должно быть совпадение. Сейчас нет совпадения ни при том а, ни при другом. То, что у тебя по -разному задается x и f, дает возможность визуально менять ширину выводимого спектра. Теоретически можно подогнать так, что разницы не будет, но вычитание векторов показывает, что они разные.
plot(f,sqrt(Puu)/max(sqrt(Puu-g, '--b', 'linewidth', 1 );

Nastia-2008

Я что-то плохо понял, что ты написал. Я не понял, почему [math]$a=\frac{1}{2\pi}$[/math].
Вообще, в матлабе подразумевается использование преобразования вида
[math]$G(f_{x})=\int \limits_{-\infty}^{\infty}g(x)e^{-i2\pi f_{x} x}dx$[/math]
и для обратного:
[math]$g(x) = \int \limits_{-\infty}^{\infty} G(f_x) e^{ i 2\pi f_x x} d f_x$[/math]
Посему функция вида [math]$e^{-ax^2}$[/math] имеет Фурье-образ [math]$\sqrt{\frac{\pi}{a}} e^{-\frac{(\pi f_x)^2}{a}}$[/math]. Отсюда функция совпадает со спектром в том случае, если [math]$a=\pi$[/math]
То, что у тебя по -разному задается x и f, дает возможность визуально менять ширину выводимого спектра.

А как они должны по-твоему задаваться? Крокодилы отдельно, яблоки отдельно :)

marandr

А как они должны по-твоему задаваться? Крокодилы отдельно, яблоки отдельно

я не против того, что они так задаются)
Просто в итоге не получается вектор, совпадающий с исходным, хотя визуально они совпадают.
Вообще, в матлабе подразумевается использование преобразования вида...

Хм, а хэлп говорит про

Иначе задать "а" я решил ибо http://ru.wikipedia.org/wiki/%D0%9F%D1%80%D0%B5%D0%BE%D0%B1%... таблица в конце, пункт 15

Nastia-2008

Хм, а хэлп говорит про
К этим выражениям из хелпа, написанные мной выше, как раз и сводятся. Если ты представишь интеграл в виде суммы от -N/2 до N/2-1, а потом переиндексируешь эту сумму, поскольку матлаб отрицательные индексы не понимает, то ты их и получишь. Только для этого надо свопнуть значения функции, что делается fftshift'ом. Об этом подробно написано, например, в книге Jason D. Schmidt "Numerical Simulation of Optical Wave Propagation With examples in MATLAB ®". Погугли её.
 
Просто в итоге не получается вектор, совпадающий с исходным, хотя визуально они совпадают.

 
plot(f,sqrt(Puu)/max(sqrt(Puu-g, '--b', 'linewidth', 1 );

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

marandr

Я понимаю, что чем шире функция, тем уже спектр. Но ты же сам пишешь, что при некотором "а" функция с ним совпадает.
Моя задача сейчас как раз в том, чтобы проверить это, понять, что алгоритм считает верно и применить его для другого импульса. Но не совпадает же...
P.S. Да, понимаю, что разные размерности, я про численное совпадение. Это лишь для проверки математики алгоритма.

Nastia-2008

Ок. Просто поверь (если не хочешь, читай книгу, которую я указывал выше что в матлабе дискретное ФП подразумевает использование формул преобразования, которые я писал в посте выше. Повторюсь.
[math]$G(f_x) = \int \limits_{-\infty}^{\infty} g(x) e^{-i 2 \pi f_x x} dx$[/math]
Допустим, ты хочешь задать ФП по другому, как указано в вики (кстати, читай английскую - там подробнее т.е.
[math]$G_2(\omega) = \frac{1}{\sqrt{2\pi}} \int \limits_{-\infty}^{\infty} g(x) e^{-i \omega x} dx$[/math]
Вспоминаем, что угловая частота [math]$\omega$[/math] и частота [math]$f_x$[/math] связаны выражением: [math]$\omega = 2 \pi f_x$[/math] Нехитрыми умозаключениями приходим к:
[math]$G_2(\omega) = \frac{1}{\sqrt{2\pi}} G(f_x)$[/math]
И, да! Таки функция [math]$e^{-x^2/2}$[/math] переходит в [math]$$e^{-\omega^2/2}$$[/math]
Как это выглядит в матлабе:
clear; clc;
L = 5; % numerical window size
N = 32; % number of samples
delta = L / N; % spatial step
delta_f = 1/(N*delta);
%numerical window in spatial and frequency domain
x = (-N/2 : N/2) * delta;
f = (-N/2 : N/2) * delta_f;
% angular frequency
omega = 2*pi*f;
%Gauss function
g2 = exp(-x.^2/2);
G2 = fftshift(fft(fftshift(g2 * delta;
figure(2)
plot(x,g2, '-g'); hold on;
plot(omega, abs(G2/sqrt(2*pi'-r')

marandr

Вот это мне очень нравится! Похоже неверно понял научника. Я подумал, что должны совпадать вектора. А видимо имелось ввиду совпадение при x=f g2=G(g2) . Попробовал на гауссе, загружаемом из файла - работает! Огромное СПАСИБО, помог разобраться)

Lene81

У тебя пространственный шаг неверно посчитан. Для N сэмплов на отрезке [0; L] delta = L/(N-1)

Nastia-2008

Да, второпях налажал.
Оставить комментарий
Имя или ник:
Комментарий: