Обратное дискретное преобразование Лапласа

@alexmonik · Регистрация: 05.12.2015

Студворк — интернет-сервис помощи студентам

Здравствуйте, уважаемые пользователи форума.
Долго думал в какой из разделов написать, так как в большей степени вопрос математики, но т.к. реалезуется всё в матлабе, поэтому спрошу здесь.

Сама проблема.
Было получено изображение в области изображений Лапаса вида:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?H(p) = \frac{p}{L{p}^{2}+pR+S}$ ,
L, R, S - действительные коэффициенты,
которое позже будет фигурировать в изображении от двух переменных. Но сначала нужно разобраться с ним.
Из данного изображения необходимо получить оригинал функции(разумеется, в дискретной форме; получать оригинал в аналитическом виде нет необходимости, да и в дальнейшем будет не целесообразно).
Т.к., в MATLAB нет готового обратного дискретного преобразования Лапласа, я так понимаю, необходимо произвести замену
$https://www.cyberforum.ru/cgi-bin/latex.cgi?p = j\omega$ ,
где j - мнимая единица, $https://www.cyberforum.ru/cgi-bin/latex.cgi?\omega$ - частота; и после этого воспользоваться обратным быстрым преобразоваением Фурье ifft(), чтобы получить оригинал изображения.
Разумеется, выполнив в лоб замену и применяя ifft(), получается ерунда.
Тренировался на злощастной синусоиде, но и с ней ничего не вышло. Далее привожу алгоритм, код и результаты.
1. взял изображение sin(x)
$https://www.cyberforum.ru/cgi-bin/latex.cgi?L{[sin(X)]} = \frac{1}{p^2+1}$

2. далее строю АЧХ (для проверки тот ли спектр получил)

Matlab M
1
2
3
4
5
6
7
f = 500; %   частота сигнала
j = sqrt(-1);% мнимая единица
Fd = 2048;  % частота дискретизации
FftL = Fd*4;    % количество линий спектра
omeg = -Fd/2:Fd/FftL:Fd/2-Fd/FftL; % частота
H = f./((j*omeg).^2+f^2); % изображение Лапласа, с заменой p = j*omeg
plot(omeg,abs(H)/2)

Обратное дискретное преобразование Лапласа

вроде всё верно

3. далее для ifft нужно произвести манипуляции со спектром, чтобы нулевая частота находилась в начале вектора частот

Matlab M
1
2
3
h1(1:FftL/2+1) = H(FftL/2:FftL); % манипуляции с отзеркаливанием спектра
h1(FftL/2+1:FftL) = H(1:FftL/2); %  ---//---
s = ifft(H,FftL);   % обратное преобразование

и вот здесь начинаются проблемы, а именно все значения s = NaN, т.е, где-то, как я понимаю, выскакиевает деление на ноль.

Буду благодарен, если ткенёте меня носом в ошибку.

и еще раз весь скрипт целиком:

Matlab M
1
2
3
4
5
6
7
8
9
10
f = 500; %   частота сигнала
j = sqrt(-1);% мнимая единица
Fd = 2048;  % частота дискретизации
FftL = Fd*4;    % количество линий спектра
omeg = -Fd/2:Fd/FftL:Fd/2-Fd/FftL; % частота
H = f./((j*omeg).^2+f^2); % изображение Лапласа, с заменой p = j*omeg
h1(1:FftL/2+1) = H(FftL/2:FftL); % манипуляции с отзеркаливанием спектра
h1(FftL/2+1:FftL) = H(1:FftL/2); %  ---//---
s = ifft(H,FftL);   % обратное преобразование
plot(omeg,abs(H)/2)

@alexmonik · 28.02.2018, 18:38 **[ТС]**

прошу прощения, там в коде опечатка, должно быть

Matlab M
1
2
s = ifft(h1,FftL);   % обратное преобразование
plot(s)

не исправил когда спектр строил.

@Annza · 01.03.2018, 13:28

Есть тонкости с этим отзеркаливанием спектра. Они связаны с тем, что ДПФ в Матлабе для четного количества отсчетов, например N = 8, выдаст вектор. 1ая компонента вектора - постоянная составляющая сигнала, со 2ой по N/2 компоненту - часть спектра в положительной области частот, N/2+1 - самосопряженный отсчет(всегда действительный) на частотах +-Fs/2, с N/2+2 по N компоненту - симметричная отрицательная часть спектра. Вот тут пытались разобраться

В вашем случае в h1 получается Inf - потому что в знаменателе L при некоторых значениях 0. Это так действительно по теории? Я просто не знаю

В общем я тут по-другому отзеркаливаю и Inf ограничиваю значением 1000:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
clear all
close all
f = 500; %   частота сигнала
j = sqrt(-1);% мнимая единица
Fd = 2048;  % частота дискретизации
FftL = Fd*4;    % количество линий спектра
omeg = -Fd/2:Fd/FftL:Fd/2-Fd/FftL; % частота
H = f./((j*omeg).^2+f^2); % изображение Лапласа, с заменой p = j*omeg
 
H = f./((j*omeg).^2+f^2); % изображение Лапласа, с заменой p = j*omeg
figure, plot(abs(H)/2)
 
h1(1:FftL/2) = H(FftL/2+1:FftL); % манипуляции с отзеркаливанием спектра
h1(FftL/2+1:FftL) = H(1:FftL/2); %  ---//---
figure, plot(abs(h1)/2)
h1(h1==Inf)=1000;
 
s = ifft(h1,FftL);   % обратное преобразование
figure, plot(s)

@alexmonik · 15.05.2018, 11:05 **[ТС]**

Итак, спустя почти 3 месяца

я решил свою проблему. Поэтому, если у кого-то есть аналогичная проблема, а именно получение из изображения в области Лапласа оригинала с помощью обратного дискретного преобразования Фурье, то вот как я это делал. (это сообщение для таких же глупеньких как я, так что за пояснение совсем элементарных вещей не пинайте)
к примеру, у нас есть сигнал вида $https://www.cyberforum.ru/cgi-bin/latex.cgi?cos(\omega t)$ длительностью [0;T] и у нас есть его изображение в области Лапласа. Теперь, из этого изображения мы хотим получить обратно косинус, но не через обратное преобразование Лапласа, а через обратное дискретное преобразования Фурье.
Первое, что стоит отметить, что табличные изображения даются на случай сигнала бесконечной длительности, поэтому:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\int_{0}^{\infty }cos(ft)exp(-st)dt = \frac{s}{s^2+f^2}$ ,
а в случае конечной длительности уже:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\int_{0}^{T }cos(ft)exp(-st)dt = \frac{s}{s^2+f^2} + \frac{fsin(fT)-scos(fT)}{s^2+f^2}exp(-sT)$
Таким образом мы получили изображение нашего сигнала. Теперь чтобы можно было применить ОДПФ, перейдём от изображения Лапласа к фурье-образу. Для чётных функций справедлива замена: $https://www.cyberforum.ru/cgi-bin/latex.cgi?s = 2\pi j\omega$ .
омега, которая пришла из подстановки $https://www.cyberforum.ru/cgi-bin/latex.cgi?s = 2\pi j\omega$ будет переменной величиной. С ней-то у меня и было больше всего проблем. Строго говоря, конечный сигнал имеет бесконечный спектр, но т.к. у нас нет бесконечного количества времени, частоту нужно брать:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\omega = [0:\frac{1}{T}:kf+1], k\geq 2$
(почему так, можно посмотреть, к примеру, в [Макс Ж. - Методы и техника обработки сигналов при физических измерениях] и [Сергиенко А.Б. - Цифровая обработка сигнала] )
И, собственно всё, можно напрямую подставлять наше изображение в функцию ifft(). Ниже код.

Matlab M
1
2
3
4
5
6
7
8
9
10
11
j1 = sqrt(-1); % мнимая единица
T = 6*pi; % длительность сигнала
f = 1; % частота сигнала
omeg = 0:1/T:20*f+1; % вектор частот
p = 2*pi*j1*omeg; % замена переменной p на j*omega
H = (p./(p.^2+f^2)) + (exp(-p*T).*(f.*sin(f*T)-p.*cos(f*T)))./(p.^2+f^2); % изображение сигнала
H(~isfinite(H))= 1; % может случиться, что числитель в верхней строке будет равен нулю
t = 0:T/(length(omeg)-1):T; % вектор времен для построения оригинала изображения
plot(t,real(ifft(H*length(omeg)))) % во втором аргументе функции plot() сразу применяется функция обратного преобразования 
                                                 %фурье ifft() и массив с изображением домножается на собственную длину (т.к. в ifft() 
                                                 %производится деление на эту самую длину)

Новые блоги и статьи Все статьи Все блоги /
Семь CDC на одном интерфейсе: 5 U[S]ARTов, 1 CAN и 1 SSI Eddy_Em 18.02.2026 Постепенно допиливаю свою "многоинтерфейсную плату". Выглядит вот так: https:/ / www. cyberforum. ru/ blog_attachment. php?attachmentid=11617&stc=1&d=1771445347 Основана на STM32F303RBT6. На борту пять. . .	Символьное дифференцирование igorrr37 13.02.2026 / * Программа принимает математическое выражение в виде строки и выдаёт его производную в виде строки и вычисляет значение производной при заданном х Логарифм записывается как: (x-2)log(x^2+2) -. . .	Камера Toupcam IUA500KMA Eddy_Em 12.02.2026 Т. к. у всяких "хикроботов" слишком уж мелкий пиксель, для подсмотра в ESPriF они вообще плохо годятся: уже 14 величину можно рассмотреть еле-еле лишь на экспозициях под 3 секунды (а то и больше),. . .	И ясному Солнцу zbw 12.02.2026 И ясному Солнцу, и светлой Луне. В мире покоя нет и люди не могут жить в тишине. А жить им немного лет.
«Знание-Сила» zbw 12.02.2026 «Знание-Сила» «Время-Деньги» «Деньги -Пуля»	SDL3 для Web (WebAssembly): Подключение Box2D v3, физика и отрисовка коллайдеров 8Observer8 12.02.2026 Содержание блога Box2D - это библиотека для 2D физики для анимаций и игр. С её помощью можно определять были ли коллизии между конкретными объектами и вызывать обработчики событий столкновения. . . .	SDL3 для Web (WebAssembly): Загрузка PNG с прозрачным фоном с помощью SDL_LoadPNG (без SDL3_image) 8Observer8 11.02.2026 Содержание блога Библиотека SDL3 содержит встроенные инструменты для базовой работы с изображениями - без использования библиотеки SDL3_image. Пошагово создадим проект для загрузки изображения. . .	SDL3 для Web (WebAssembly): Загрузка PNG с прозрачным фоном с помощью SDL3_image 8Observer8 10.02.2026 Содержание блога Библиотека SDL3_image содержит инструменты для расширенной работы с изображениями. Пошагово создадим проект для загрузки изображения формата PNG с альфа-каналом (с прозрачным. . .

Обратное дискретное преобразование Лапласа

Решение