diff --git a/.github/workflows/build-pdf.yml b/.github/workflows/build-pdf.yml index d0d7f31..fa6a5ab 100644 --- a/.github/workflows/build-pdf.yml +++ b/.github/workflows/build-pdf.yml @@ -29,12 +29,16 @@ jobs: 3_fft_selfmade.md \ 4_ifft_selfmade.md \ 5_comparison.md \ + 6_additional_exs.md \ -o coursework.pdf \ --pdf-engine=xelatex \ -V mainfont="Liberation Serif" \ + -V monofont="Liberation Mono" \ -V lang=ru-RU \ -V geometry:margin=2cm \ - --toc + --toc \ + --toc-depth=2 \ + --number-sections - name: Upload PDF uses: actions/upload-artifact@v4 diff --git a/.github/workflows/build-presentation.yml b/.github/workflows/build-presentation.yml new file mode 100644 index 0000000..0d36cef --- /dev/null +++ b/.github/workflows/build-presentation.yml @@ -0,0 +1,35 @@ +name: Build Marp Presentations (PDF) + +on: + push: + branches: [ "main", "master" ] + paths: + - 'slides/**.md' + workflow_dispatch: + +jobs: + build: + runs-on: ubuntu-latest + + steps: + - name: checkout + uses: actions/checkout@v4 + + - name: install node + uses: actions/setup-node@v4 + with: + node-version: '20' + + - name: generate with marp + run: | + npx @marp-team/marp-cli@latest slides/ --pdf --allow-local-files + env: + CHROME_PATH: "" + PUPPETEER_SKIP_CHROMIUM_DOWNLOAD: "false" + + - name: save artifacts + uses: actions/upload-artifact@v4 + with: + name: Coursework-Presentations + path: slides/*.pdf + retention-days: 30 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..cf9f386 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +venv +octave-workspace +.DS_Store \ No newline at end of file diff --git a/.zprofile b/.zprofile index 0e9cafc..9d94a73 100644 --- a/.zprofile +++ b/.zprofile @@ -1,3 +1,108 @@ export APT_HOME=~/Downloads/APT_v3.1.0 export PATH=$PATH:$APT_HOME export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC +export APT_HOME=~/Downloads/APT_v3.1.0 +export PATH=$PATH:$APT_HOME +export APT_ARCH=MAC diff --git a/1_analytical.md b/1_analytical.md new file mode 100644 index 0000000..6db9c1a --- /dev/null +++ b/1_analytical.md @@ -0,0 +1,133 @@ +# Самостоятельная работа: Дискретное преобразование Фурье. Вариант 1 + +## Аналитический расчет коэффициентов ряда Фурье + + +$$ +\left\{ +\begin{aligned} +a_0 &= 0.1 \text{ МПа} \\ +f_0 &= 2 \text{ МГц} \\ +\omega_0 &= 2\pi f_0 \\ +p(t) &= 2a_0 \sin(3\omega_0 t) + a_0 \cos(5\omega_0 t) +\end{aligned} +\right. +$$ + +### Определение периода $T$ +Сигнал состоит из двух гармоник с частотами $f_1 = 3f_0 = 6$ МГц и $f_2 = 5f_0 = 10$ МГц. + +Определим базовую частоту: +$$ +f_0 = \text{НОД(6 МГц, 10 МГц)} = \text{2МГц} +$$ + + +Следовательно, период сигнала $T$ равен: +$$ +T = \frac{1}{f_0} = \frac{1}{2 \cdot 10^6} = 0.5 \text{ мкс} +$$ + +### Расчет аналитических коэффициентов спектра $p_T(f_n)$ +Используем формулы Эйлера для представления тригонометрических функций через комплексные экспоненты: + +$$ +p(t) = 2a_0 \sin(3\omega_0 t) + a_0 \cos(5\omega_0 t) +$$ + + +$$ +\implies +\left[ +\begin{aligned} + \sin(x) = \frac{e^{ix} - e^{-ix}}{2i} = -\frac{i}{2}e^{ix} + \frac{i}{2}e^{-ix} \\ + \cos(x) = \frac{e^{ix} + e^{-ix}}{2} = \frac{1}{2}e^{ix} + \frac{1}{2}e^{-ix} +\end{aligned} +\right] +$$ +$$ \implies p(t) = 2(0.1) \left[ -\frac{i}{2}e^{i 3\omega_0 t} + \frac{i}{2}e^{-i 3\omega_0 t} \right] + 0.1 \left[ \frac{1}{2}e^{i 5\omega_0 t} + \frac{1}{2}e^{-i 5\omega_0 t} \right] $$ +$$ \implies p(t) = -0.1i \cdot e^{i 3\omega_0 t} + 0.1i \cdot e^{-i 3\omega_0 t} + 0.05 \cdot e^{i 5\omega_0 t} + 0.05 \cdot e^{-i 5\omega_0 t} $$ + +Из полученного выражения напрямую следуют аналитические значения комплексных коэффициентов спектра $p_T(f_n)$ на частотах $f_n$: + + +$$ +\begin{aligned} +& 3f_0 (6 МГц): p_T(3f_0) = -0.1i & \text{Фаза}& = -\pi/2 \\ +& -3f_0 (-6 МГц): p_T(-3f_0) = 0.1i; & \text{Фаза}&= \pi/2 \\ +&5f_0 (10 МГц): p_T(5f_0) = 0.05; & \text{Фаза}& = 0 \\ +&-5f_0 (-10 МГц): p_T(-5f_0) = 0.05; & \text{Фаза}&= 0\\ +\end{aligned} +$$ + +### Графики функции и аналитического спектра + +Амплитудный и действительный спектры четные (рис. 1.a, 1.c), а фазовый и мнимый спектры нечетные (рис. 1.b, 1.d). + +![Аналитический спектр и сигнал](pic/fig_1_m.png) + +\newpage + +### Код №1 +```matlab +a_0 = 0.1; +f_0 = 2.0; +w_0 = 2 * 3.14 * f_0; +T = 1 / f_0; +t = 0: 2*T: 1000; +p_t = 2 * a_0 * sin(3*w_0 * t) + a_0 * cos(5 * w_0 * t); +freqs_f0_multipliers = -6: 1: 6; +freqs = freqs_f0_multipliers * f_0; + +coefficients = zeros(length(freqs), length(t)); +coefficients(freqs_f0_multipliers == 3) = -0.1i; +coefficients(freqs_f0_multipliers == -3) = 0.1i; +coefficients(freqs_f0_multipliers == 5) = 0.05; +coefficients(freqs_f0_multipliers == -5) = 0.05; + +c_amp = abs(coefficients); +c_phase = angle(coefficients); +c_real = real(coefficients); +c_im = imag(coefficients); + +figure('Position', [110, 55, 1200, 1200]); + +axes('Position', [0.0667, 0.7066, 0.9208, 0.2606]); +plot(t, p_t); +title('Исходный акустический сигнал p(t)'); +xlabel('Время t, мкс'); +ylabel('p(t), МПа'); +grid on; + +axes('Position', [0.08, 0.38, 0.39, 0.22]); +stem(freqs, c_amp); +hold on; plot(xlim, [0 0]); hold off; +title('Амплитудный спектр |p_T(f_n)|'); +xlabel('Частота f_n, МГц'); +grid on; + +axes('Position', [0.57, 0.38, 0.39, 0.22]); +stem(freqs, c_phase); +hold on; plot(xlim, [0 0]); hold off; +title('Фазовый спектр arg(p_T(f_n))') +xlabel('Частота f_n, МГц'); +ylabel('Фаза p_T(f_n)'); +grid on; + +axes('Position', [0.08, 0.07, 0.39, 0.22]); +stem(freqs, c_real); +hold on; plot(xlim, [0 0]); hold off; +title('Действительная часть спектра Re(p_T(f_n))') +xlabel('Частота f_n, МГц'); +grid on; + +axes('Position', [0.57, 0.07, 0.39, 0.22]); +stem(freqs, c_im); +hold on; plot(xlim, [0 0]); hold off; +title('Мнимая часть спектра Im(p_T(f_n))') +xlabel('Частота f_n, МГц'); +grid on; + + +print(gcf, 'fig_1_m.png', '-dpng', '-r300'); +``` diff --git a/2_discrette.md b/2_discrette.md new file mode 100644 index 0000000..6944598 --- /dev/null +++ b/2_discrette.md @@ -0,0 +1,33 @@ +## Выбор параметров дискретизации + +### Теорема Котельникова и шаг дискретизации $h$ +Согласно теореме Котельникова, частота дискретизации $f_s$ должна быть более чем в два раза выше максимальной частоты, присутствующей в спектре сигнала: +$$ f_s > 2 f_{max} $$ + +Из результатов аналитического расчета: +$$ f_{max} = 5f_0 = 10 \text{ МГц} \implies f_s > 20 \text{ МГц} $$ + +Следовательно, минимально допустимая частота дискретизации: +$$ f_{s_{min}} = 2 \cdot 10 \text{ МГц} = 20 \text{ МГц} $$ + +Для шага дискретизации $h$ выберем максимально допустимое значение: +$$ h < \frac{1}{20 \text{ МГц}} = 0.05 \text{ мкс} $$ + +### Выбор временного окна $T$ +Для ДПФ периодического сигнала размер окна должен быть кратным периоду сигнала, поэтому выберем временное окно равным ровно одному периоду сигнала: +$$ T_{window} = T = 0.5 \text{ мкс} $$ + +### Определение количества точек $N$ +Рассчитаем необходимое минимальное количество точек $N^*$ для нашего окна при максимальном шаге дискретизации: +$$ N^* = \frac{T_{window}}{h_{max}} = \frac{0.5}{0.05} = 10 \text{ точек} $$ + +Чтобы пакетная реализация fft работала эффективно, выберем количество точек как ближайшую степень двойки, превышающую $N^*$: + +$$ \mathbf{N = 16} $$ + +### Итоговые параметры расчетной сетки +Пересчитаем фактический шаг дискретизации $h$ для выбранного $N$: +$$ h = \frac{T_{window}}{N} = \frac{0.5}{16} = 0.03125 \text{ мкс} $$ + +Частота дискретизации при этом составит: +$$ f_s = \frac{1}{h} = 32 \text{ МГц} $$ diff --git a/3_fft_selfmade.md b/3_fft_selfmade.md new file mode 100644 index 0000000..763930f --- /dev/null +++ b/3_fft_selfmade.md @@ -0,0 +1,101 @@ +## Программная реализация ДПФ + +Реализуем алгоритм ДПФ для вычисления спектра сеточной функции $p(l)$ по формуле из задания: + +$$ +\overline{p_T(n)} = \frac{1}{N} \sum_{l=0}^{N-1} p(l)W_N^{-nl}, \quad W_N^{nl} = \exp\left( i\frac{2\pi nl}{N} \right) \qquad +$$ + + +### Использование свойств симметрии +Поскольку исходный сигнал $p(t)$ является действительной функцией, его спектр обладает свойством эрмитовой симметрией. Для $N$-точечного ДПФ: +$$ p_T(N - n) = p_T^*(n) $$ +где $*$ обозначает комплексное сопряжение. + +Будем сначала вычислять спектральные коэффициенты $p_T(n)$ по формуле суммы только для первой половины спектра (от $n = 0$ до $n = N/2$), а вторую половину получим операцией комплексного сопряжения. + + +### Построение графиков +Графики построены в расширенном интервале от $-1/h$ до $1/h$ (т.е. от $-f_s$ до $f_s$). Видно, что амплитуды получились в два раза меньше рассчитанных аналитически (рис. 3.a). Это следствие того, что мы рассматриваем интервал, включающий отрицательные частоты. +Например: $$A_{analytical}⋅cos⁡(5ω_0t)=\frac{A_{analytical}}{2}⋅e^{i5ω0t}+\frac{A_{analytical}}{2}⋅e^{−i5ω0t}$$ + +![Собственное ДПФ](pic/fig_3_m.png) + +\newpage + +### Код №3 + +```matlab +a0 = 0.1; +f0 = 2.0; +w0 = 2 * pi * f0; +T = 0.5; +N = 16; +h = T / N; +fs = 1 / h; + +l_indices = 0:1:N; +t_l = l_indices * h; +p_l = 2 * a0 * sin(3*w0 * t_l) + a0 * cos(5 * w0 * t_l); +p_T_discrete = zeros(N); + +half_N = N / 2; +for n = 0 : half_N + sum_val = 0; + for l = 0:(N-1) + exponent = -1i * 2 * pi * n * l / N; + sum_val = sum_val + p_l(l+1) * exp(exponent); + end + p_T_discrete(n+1) = sum_val / N; +end + +for n = half_N+1:(N-1) + p_T_discrete(n) = conj(p_T_discrete((N - n) + 1)); +end + +n_extended = -N:1:N; +freqs_extended = n_extended .* (fs / N); +p_T_extended = p_T_discrete(mod(n_extended, N) + 1); +p_T_extended(abs(p_T_extended) < 1e-10) = 0; + +ampliude_ext = abs(p_T_extended); +phase_ext = angle(p_T_extended); +real_part_ext = real(p_T_extended); +imag_part_ext = imag(p_T_extended); + + +figure('Position', [110, 55, 1600, 1300]); + +axes('Position', [0.0628, 0.5508, 0.4298, 0.4129]); +stem(freqs_extended, ampliude_ext); +hold on; +y_bounds = ylim; +plot([fs/2, fs/2], y_bounds); +plot([-fs/2, -fs/2], y_bounds); +hold off; +title('Амплитудный спектр ДПФ'); +xlabel('Частота, МГц'); +ylabel('Амплитуда, МПа'); +grid on; + +axes('Position', [0.5577, 0.5508, 0.4298, 0.4129]); +stem(freqs_extended, phase_ext) +title('Фазовый спектр ДПФ'); +xlabel('Частота, МГц'); +ylabel('Фаза, рад'); +grid on; + +axes('Position', [0.0628, 0.0583, 0.4298, 0.4129]); +stem(freqs_extended, real_part_ext) +title('Действительная часть ДПФ'); +xlabel('Частота, МГц'); +grid on; + +axes('Position', [0.5577, 0.0583, 0.4298, 0.4129]); +stem(freqs_extended, imag_part_ext) +title('Мнимая часть ДПФ'); +xlabel('Частота, МГц'); +grid on; + +print(gcf, 'fig_3_m.png', '-dpng', '-r300'); +``` \ No newline at end of file diff --git a/4_ifft_selfmade.md b/4_ifft_selfmade.md new file mode 100644 index 0000000..e9a78ac --- /dev/null +++ b/4_ifft_selfmade.md @@ -0,0 +1,112 @@ +## Обратное ДПФ и интерполяция по теореме Котельникова + +В этом разделе решается обратная задача: восстановление временного сигнала из его дискретного спектра (ОДПФ) и реконструкция непрерывной формы сигнала между узлами дискретизации. + +### Программная реализация обратного ДПФ (ОДПФ) +Вычисление сеточной функции $p(l)$ по найденным амплитудам гармоник проводится по формуле: +$$ p(l) = \sum_{n=0}^{N-1} p_T(n) \exp\left( i \frac{2\pi nl}{N} \right) $$ + +Теоретически, прямое и обратное преобразования формируют строго обратимую пару. Это можно видеть на графиках в конце данного раздела. + +### Восстановление непрерывного сигнала + +Для восстановления непрерывной функции $p(t)$ используется ряд Котельникова (формула 5): +$$ p(t) = \sum_{l=0}^{N-1} p(l) \text{sinc} \left( \frac{\pi}{h} (t - lh) \right) $$ + +В библиотеке `NumPy` и в реализации sinc в `Matlab` используется нормированная функция `sinc(x)`, которая вычисляется как $\frac{\sin(\pi x)}{\pi x}$. Поэтому при программировании аргумент функции умножать на $\pi$ не нужно, формула принимает вид `numpy.sinc((t - l*h) / h)`. + +### Графики результатов +Вычисления производились на сетке с шагом $h/10$. Видно, что интерполированная кривая проходит через исходные узлы дискретизации и совпадает с исходной аналитической функцией, что доказывает теорему Котельникова (рис. 4.a, 4.b). + +![ОДПФ и интерполяция](pic/fig_4_m.png) + + +### Оценка ошибки + +Видно, что для каждой из точек ошибка не превышает машинную точность (порядка $10^{-16}$). + +![Оценка ошибки](pic/fig_4_error.png) + +\newpage + +### Код № 4 + +```matlab +a0 = 0.1; +f0 = 2.0; +w0 = 2 * pi * f0; +T = 0.5; +N = 16; +h = T / N; + +l_indices = 0 : (N - 1); +t_l = l_indices * h; +p_l = 2 * a0 * sin(3 * w0 * t_l) + a0 * cos(5 * w0 * t_l); + +p_T_discrete = zeros(1, N); +for n = 0 : (N - 1) + sum_val = 0; + for l = 0 : (N - 1) + exponent = -1i * 2 * pi * n * l / N; + sum_val = sum_val + p_l(l + 1) * exp(exponent); + end + p_T_discrete(n + 1) = sum_val / N; +end + +p_restored = zeros(1, N); +for l = 0 : (N - 1) + sum_val = 0; + for n = 0 : (N - 1) + exponent = 1i * 2 * pi * n * l / N; + sum_val = sum_val + p_T_discrete(n + 1) * exp(exponent); + end + p_restored(l + 1) = sum_val; +end + +p_restored_real = real(p_restored); + +h_fine = h / 10; + +t_fine = 0 : h_fine : (T - h_fine); + +p_interpolated = zeros(1, length(t_fine)); + +for idx = 1 : length(t_fine) + t = t_fine(idx); + sum_sinc = 0; + for l = 0 : (N - 1) + sum_sinc = sum_sinc + p_l(l + 1) * sinc((t - l * h) / h); + end + p_interpolated(idx) = sum_sinc; +end + +p_analytical_fine = 2 * a0 * sin(3 * w0 * t_fine) + a0 * cos(5 * w0 * t_fine); + +figure('Position', [110, 55, 1500, 1050]); + +axes('Position', [0.0781, 0.5651, 0.9069, 0.3890]); +plot(t_l, p_l, 'Исходный дискретный p(l)'); +hold on; +plot(t_l, p_restored_real, 'Восстановленный ОДПФ'); +hold off; + +title('Проверка обратимости прямого и обратного ДПФ'); +xlabel('Время t, мкс'); +ylabel('Амплитуда, МПа'); +grid on; +legend('Location', 'best'); + +axes('Position', [0.0781, 0.0747, 0.9069, 0.3890]); +plot(t_l, p_l, 'ko', 'Узлы дискретизации p(l)'); +hold on; +plot(t_fine, p_interpolated, 'Интерполяция (sinc)'); +plot(t_fine, p_analytical_fine, 'Аналитический сигнал'); +hold off; +title('Восстановление непрерывного сигнала (Теорема Котельникова)'); +xlabel('Время t, мкс'); +ylabel('Амплитуда, МПа'); +grid on; +legend('Location', 'best'); + +print(gcf, 'fig_4_m.png', '-dpng', '-r300'); +``` diff --git a/5_comparison.md b/5_comparison.md new file mode 100644 index 0000000..fbc6ddc --- /dev/null +++ b/5_comparison.md @@ -0,0 +1,76 @@ +## Использование встроенных функций БПФ (FFT) + +В данном разделе для расчета спектра и восстановления сигнала применяются встроенные библиотечные функции быстрого преобразования Фурье (`fft`) и обратного быстрого преобразования Фурье (`ifft`) из `Matlab`. + +Из-за особенности реализации `fft`, соотношение между встроенным БПФ и формулой из задания №3 имеет вид: +$$ \text{fft}(p(l)) = N \cdot p_T(n) $$ + +Чтобы получить физические амплитуды гармоник, результат функции `fft` необходимо разделить на размер выборки $N$. + +### Проверка обратимости БПФ +Для проверки обратимости выполним обратное преобразование полученного БПФ-спектра: +$$ p_{IFFT} = \text{ifft}(P_{FFT}) $$ + +Результаты вычислений и график показывают, что восстановленный сигнал $p_{IFFT}$ точно совпадает с исходной сеточной функцией $p(l)$. + +Алгоритмы `fft` и `ifft` составляют строго обратимую пару. + +![Результаты работы БПФ](pic/fig_5_m.png) + + +### Оценка ошибки + +Видно, что ошибка не превышает машинную точность (порядка $10^{-17}$), также на один порядок меньше ошибки в собственной реализации. + +![Оценка ошибки](pic/fig_5_error.png) + + +\newpage + +### Код + +```matlab +a0 = 0.1; +f0 = 2.0; +w0 = 2 * pi * f0; +T = 0.5; +N = 16; +h = T / N; +l_indices = 0 : (N - 1); +t_l = l_indices * h; +p_l = 2 * a0 * sin(3 * w0 * t_l) + a0 * cos(5 * w0 * t_l); +P_fft_raw = fft(p_l); +P_fft_normalized = P_fft_raw / N; +amplitude_fft = abs(P_fft_normalized); +fs = 1 / h; +freqs_fft = (0 : (N - 1)) * (fs / N); +p_ifft = ifft(P_fft_raw); +p_ifft_real = real(p_ifft); + + +figure('Position', [110, 55, 1600, 1300]); +axes('Position', [0.0650, 0.5513, 0.9225, 0.4114]); +half_N = (N / 2) + 1; +stem(freqs_fft(1:half_N), amplitude_fft(1:half_N), 'Спектр ДПФ'); +hold on; +plot(6.0, 0.1, 'ro', 'Аналитическая амплитуда (3f_0)'); +plot(10.0, 0.05, 'ro', 'Аналитическая амплитуда (5f_0)'); +hold off; +title('Амплитудный спектр, полученный с помощью БПФ (нормированный на N)'); +xlabel('Частота, МГц'); +ylabel('Амплитуда, МПа'); +grid on; +legend('Location', 'northeast'); + +axes('Position', [0.0650, 0.0598, 0.9225, 0.4114]); +plot(t_l, p_l, 'k-', 'LineWidth', 2, 'DisplayName', 'Исходный сигнал p(l)'); +hold on; +plot(t_l, p_ifft_real, 'Восстановленный ifft'); +hold off; +title('Проверка обратимости прямого и обратного БПФ (FFT / IFFT)'); +xlabel('Время \itt\rm, мкс'); +ylabel('Амплитуда, МПа'); +grid on; +legend('Location', 'northeast'); +print(gcf, 'fig_5_m.png', '-dpng', '-r300'); +``` diff --git a/6_additional_exs.md b/6_additional_exs.md new file mode 100644 index 0000000..b45b077 --- /dev/null +++ b/6_additional_exs.md @@ -0,0 +1,162 @@ +# Дополнительное задание + +## Аналитический расчет спектра производной + +По свойствам преобразования Фурье, дифференцирование сигнала во временной области эквивалентно умножению его спектра на $i \omega$ (или $i 2\pi f_n$) в частотной области: +$$ \frac{dp(t)}{dt} \longleftrightarrow i 2\pi f_n \cdot P(f_n) = i \omega_n \cdot P(\omega_n) $$ + +Найдем производную аналитически: +$$ +p(t) = 2a_0 \sin(3\omega_0 t) + a_0 \cos(5\omega_0 t) = 0.2 \sin(3\omega_0 t) + 0.1 \cos(5\omega_0 t) +$$ + +$$ \implies \frac{dp(t)}{dt} = 0.2 \cdot 3\omega_0 \cos(3\omega_0 t) - 0.1 \cdot 5\omega_0 \sin(5\omega_0 t) = 0.6\omega_0 \cos(3\omega_0 t) - 0.5\omega_0 \sin(5\omega_0 t) $$ + +Применим свойство к аналитическим коэффициентам из п.1.: + +$$ +\begin{aligned} +& 3f_0 (6 МГц): C'_3 = i(2\pi \cdot 3f_0) \cdot (-0.1i) = 0.6\pi f_0 = 1.2\pi \approx 3.77 & \text{(Чисто действительное)} \\ +& -3f_0 (-6 МГц): C'_{-3} = i(2\pi \cdot (-3f_0)) \cdot (0.1i) = 0.6\pi f_0 = 1.2\pi \approx 3.77 \\ +&5f_0 (10 МГц): C'_5 = i(2\pi \cdot 5f_0) \cdot 0.05 = i \cdot 0.5\pi f_0 = i \pi \approx 3.14 i & \text{(Чисто мнимое)} \\ +&-5f_0 (-10 МГц): C'_{-5} = i(2\pi \cdot (-5f_0)) \cdot 0.05 = -i \cdot 0.5\pi f_0 = -i \pi \approx -3.14 i\\ +\end{aligned} +$$ + + +Производная действительного сигнала также является действительным сигналом. Cпектр обладает свойством эрмитовой симметрии: действительная часть спектра и амплитуда — четные функции ($C'_3 = C'_{-3}$), а мнимая часть и фаза — нечетные функции ($C'_5 = -C'_{-5}$). Умножение на $i$ привело к тому, что гармоники, бывшие чисто мнимыми, стали действительными, и наоборот (рис. 6). + +![Спектр производной (аналитический)](pic/fig_add_1_m.png) + +## Расчет дискретного спектра производной из спектра сигнала + +Для получения спектра производной из дискретного спектра $P_{discrete}(n)$ (рассчитанного в п. 3) необходимо каждый отсчет умножить на $i 2\pi f_n$. + +В дискретном спектре индексы $n$ от $0$ до $N/2$ соответствуют положительным частотам, а индексы от $N/2$ до $N-1$ представляют область отрицательных частот. Чтобы корректно применить множитель $i \omega$, необходимо ввести масштабирование частоты — пересчитать индекс $n$ в физическую частоту: +* Для $0 \le n \le N/2$: $f_{mapped} = n \cdot \frac{f_s}{N}$ +* Для $N/2 < n \le N-1$: $f_{mapped} = (n - N) \cdot \frac{f_s}{N}$ + +Таким образом, спектр производной рассчитывается как: +$$ P'_{discrete}(n) = P_{discrete}(n) \cdot \left( i 2\pi f_{mapped}(n) \right) $$ + +Видно, что амплитуды на графиках совпадают с п.1, что подтверждает правильность введенного частотного масштабирования. + +![Дискретный спектр производной](pic/fig_add_2_m.png) + +## Расчет спектра производной с использованием БПФ (FFT / IFFT) + +В данном разделе применяется описанный выше алгоритм, но с использованием встроенных функций MATLAB/Python: +1. Вычисляется спектр сигнала: $P = \text{fft}(p)$ +2. Формируется вектор круговых частот $\omega_{vec}$ с учетом отрицательной полуоси (функция `fftfreq`). +3. Вычисляется спектр производной: $P'_{FFT} = P \cdot (i \omega_{vec})$ +4. Восстанавливается сигнал во временной области: $p'_{IFFT} = \text{ifft}(P'_{FFT})$ + +На графике ниже представлено сравнение восстановленной через IFFT производной и точного аналитического решения для $dp/dt$. + + +![БПФ спектральное дифференцирование](pic/fig_add_3_m.png) + + +Также стоит заметить, что ошибка вычисления не превышает машинную точность вычислений (порядка $10^{-14}$). + +![Оценка ошибки](pic/fig_add_6_error.png) + +\newpage + +### Код + +```matlab +a0 = 0.1; +f0 = 2.0; +w0 = 2 * pi * f0; +T = 0.5; +N = 16; +h = T / N; +fs = 1 / h; +l_indices = 0 : (N - 1); +t_l = l_indices * h; +p_l = 2 * a0 * sin(3 * w0 * t_l) + a0 * cos(5 * w0 * t_l); +freqs_f0_multipliers = -6 : 6; +freqs = freqs_f0_multipliers * f0; +coeffs_deriv = zeros(1, length(freqs)); +coeffs_deriv(freqs_f0_multipliers == 3) = 0.6 * pi * f0; +coeffs_deriv(freqs_f0_multipliers == -3) = 0.6 * pi * f0; +coeffs_deriv(freqs_f0_multipliers == 5) = 1i * 0.5 * pi * f0; +coeffs_deriv(freqs_f0_multipliers == -5) = -1i * 0.5 * pi * f0; +figure('Position', [110, 55, 1600, 1300]); +axes('Position', [0.0391, 0.5390, 0.4498, 0.4156]); +stem(freqs, abs(coeffs_deriv), 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +title('Амплитудный спектр dp/dt'); +axes('Position', [0.5377, 0.5390, 0.4498, 0.4156]); +stem(freqs, angle(coeffs_deriv), 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +title('Фазовый спектр dp/dt'); +axes('Position', [0.0391, 0.0484, 0.4498, 0.4156]); +stem(freqs, real(coeffs_deriv), 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +title('Действительная часть dp/dt'); +axes('Position', [0.5377, 0.0484, 0.4498, 0.4156]); +stem(freqs, imag(coeffs_deriv), 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +title('Мнимая часть dp/dt'); +print(gcf, 'fig_add_1_m.png', '-dpng', '-r300'); +p_T_discrete = zeros(1, N); +for n = 0 : (N - 1) + sum_val = 0; + for l = 0 : (N - 1) + exponent = -1i * 2 * pi * n * l / N; + sum_val = sum_val + p_l(l + 1) * exp(exponent); + end + p_T_discrete(n + 1) = sum_val / N; +end +freqs_mapped = zeros(1, N); +for n = 0 : (N - 1) + if n <= N/2 + freqs_mapped(n + 1) = n * (fs / N); + else + freqs_mapped(n + 1) = (n - N) * (fs / N); + end +end +omega_mapped = 2 * pi * freqs_mapped; +p_T_deriv_discrete = p_T_discrete .* (1i * omega_mapped); +n_extended = -N : N; +freqs_extended = n_extended * (fs / N); +p_T_deriv_extended = p_T_deriv_discrete(mod(n_extended, N) + 1); +p_T_deriv_extended(abs(p_T_deriv_extended) < 1e-10) = 0; + +figure('Position', [110, 55, 1500, 1050]); +axes('Position', [0.0391, 0.5390, 0.4498, 0.4141]); +stem(freqs_extended, abs(p_T_deriv_extended)); +hold on; +y_bounds = ylim; +plot([fs/2, fs/2], y_bounds); +plot([-fs/2, -fs/2], y_bounds); +hold off; +title('Амплитуда P_discrete в интервале [-1/h, 1/h]'); +axes('Position', [0.5377, 0.5390, 0.4498, 0.4141]); +stem(freqs_extended, angle(p_T_deriv_extended)); +title('Фаза P_{discrete}'); +axes('Position', [0.0391, 0.0484, 0.4498, 0.4141]); +stem(freqs_extended, real(p_T_deriv_extended)); +title('Re(P_{discrete})'); +axes('Position', [0.5377, 0.0484, 0.4498, 0.4141]); +stem(freqs_extended, imag(p_T_deriv_extended)); +title('Im(P_{discrete})'); +print(gcf, 'fig_add_2_m.png', '-dpng', '-r300'); +P_fft = fft(p_l); +% there is no fftfreq in matlab, +% so these freqs came out from python script +freqs_fft = [0 : (N/2 - 1), -N/2 : -1] * (fs / N); +omega_fft = 2 * pi * freqs_fft; +P_deriv_fft = P_fft .* (1i * omega_fft); +p_deriv_ifft = ifft(P_deriv_fft); +p_deriv_ifft_real = real(p_deriv_ifft); +p_deriv_analytical = 0.6 * w0 * cos(3 * w0 * t_l) - 0.5 * w0 * sin(5 * w0 * t_l); +figure('Position', [110, 55, 1500, 750]); +plot(t_l, p_deriv_analytical, 'Аналитическая dp/dt'); +hold on; +plot(t_l, p_deriv_ifft_real, 'Восстановленная ifft'); +hold off; +title('Восстановление производной через обратное БПФ'); +xlabel('Время t, мкс'); +ylabel('Амплитуда, МПа/мкс'); +legend('Location', 'northeast'); +print(gcf, 'fig_add_3_m.png', '-dpng', '-r300'); +``` diff --git a/code/1.m b/code/1.m new file mode 100644 index 0000000..b2313ac --- /dev/null +++ b/code/1.m @@ -0,0 +1,85 @@ +clear; clc; close all; + +blue_color = [0.12, 0.47, 0.71]; +red_color = [0.84, 0.15, 0.16]; +set(0, 'DefaultFigureColor', 'w'); +set(0, 'DefaultAxesColor', 'w'); +set(0, 'DefaultTextFontName', 'DejaVu Sans'); +set(0, 'DefaultAxesFontName', 'DejaVu Sans'); +set(0, 'DefaultAxesFontSize', 10); +set(0, 'DefaultTextFontSize', 12); +set(0, 'DefaultTextFontWeight', 'normal'); +set(0, 'DefaultAxesTitleFontWeight', 'normal'); +set(0, 'DefaultAxesTickDir', 'in'); +set(0, 'DefaultAxesLineWidth', 1.5); +set(0, 'DefaultLineLineWidth', 1.7); +set(0, 'DefaultAxesTickLength', [0.015 0.025]); + +set(0, 'DefaultAxesXGrid', 'on'); +set(0, 'DefaultAxesYGrid', 'on'); +set(0, 'DefaultAxesGridLineStyle', '-'); +set(0, 'DefaultAxesGridColor', [0.6 0.6 0.6]); +set(0, 'DefaultAxesGridAlpha', 0.5); + +a_0 = 0.1; +f_0 = 2.0; +w_0 = 2 * 3.14 * f_0; +T = 1 / f_0; +t = 0: 2*T: 1000; +p_t = 2 * a_0 * sin(3*w_0 * t) + a_0 * cos(5 * w_0 * t); +freqs_f0_multipliers = -6: 1: 6; +freqs = freqs_f0_multipliers * f_0; + +coefficients = zeros(length(freqs), length(t)); +coefficients(freqs_f0_multipliers == 3) = -0.1i; +coefficients(freqs_f0_multipliers == -3) = 0.1i; +coefficients(freqs_f0_multipliers == 5) = 0.05; +coefficients(freqs_f0_multipliers == -5) = 0.05; + +c_amp = abs(coefficients); +c_phase = angle(coefficients); +c_real = real(coefficients); +c_im = imag(coefficients); + +c_amp(c_amp == 0) = NaN; +c_phase(c_phase == 0) = NaN; +c_real(c_real == 0) = NaN; +c_im(c_im == 0) = NaN; + +figure('Position', [110, 55, 1200, 1200]); + +axes('Position', [0.0667, 0.7066, 0.9208, 0.2606]); +plot(t, p_t, 'b', 'LineWidth', 1.7); +title('Исходный акустический сигнал p(t)'); +xlabel('Время t, миллисекунды'); +ylabel('p(t), МПа'); +grid on; + + +axes('Position', [0.08, 0.38, 0.39, 0.22]); +stem(freqs, c_amp, 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +title('\bf a) \rm Амплитудный спектр |p_T(f_n)|'); +xlabel('Частота f_n, МГц'); +grid on; + +axes('Position', [0.57, 0.38, 0.39, 0.22]); +stem(freqs, c_phase, 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +title('\bf b) \rm Фазовый спектр arg(p_T(f_n))') +xlabel('Частота f_n, МГц'); +ylabel('Фаза p_T(f_n)'); +grid on; + +axes('Position', [0.08, 0.07, 0.39, 0.22]); +stem(freqs, c_real, 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +title('\bf c) \rm Действительная часть спектра Re(p_T(f_n))') +xlabel('Частота f_n, МГц'); +grid on; + +axes('Position', [0.57, 0.07, 0.39, 0.22]); +stem(freqs, c_im, 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color) +title('\bf d) \rm Мнимая часть спектра Im(p_T(f_n))') +xlabel('Частота f_n, МГц'); +grid on; + + +print(gcf, 'fig_1_m.png', '-dpng', '-r300'); diff --git a/code/1.py b/code/1.py new file mode 100644 index 0000000..7e12fba --- /dev/null +++ b/code/1.py @@ -0,0 +1,67 @@ +import numpy as np +import matplotlib.pyplot as plt + +a0 = 0.1 # MPa +f0 = 2.0 # MHz +w0 = 2 * np.pi * f0 +T = 1 / f0 # microseconds + +t = np.linspace(0, 2*T, 1000) +p_t = 2 * a0 * np.sin(3 * w0 * t) + a0 * np.cos(5 * w0 * t) + +freqs_f0_multipliers = np.arange(-6, 7) +freqs = freqs_f0_multipliers * f0 + +coeffs = np.zeros(len(freqs), dtype=complex) + +coeffs[freqs_f0_multipliers == 3] = -0.1j +coeffs[freqs_f0_multipliers == -3] = 0.1j +coeffs[freqs_f0_multipliers == 5] = 0.05 +coeffs[freqs_f0_multipliers == -5] = 0.05 + +amplitude = np.abs(coeffs) +phase = np.angle(coeffs) +real_part = np.real(coeffs) +imag_part = np.imag(coeffs) + +fig = plt.figure(figsize=(12, 12)) + +plt.subplot(3, 1, 1) +plt.plot(t, p_t, 'b-', linewidth=2) +plt.title('Исходный акустический сигнал $p(t)$', fontsize=14) +plt.xlabel('Время $t$, мкс') +plt.ylabel('$p(t)$, МПа') +plt.grid(True) + + +plt.subplot(3, 2, 3) +plt.stem(freqs, amplitude) +plt.title('Амплитудный спектр $|p_T(f_n)|$') +plt.xlabel('Частота, МГц') +plt.grid(True) + +plt.subplot(3, 2, 4) +plt.stem(freqs, phase) +plt.title('Фазовый спектр $arg(p_T(f_n))$') +plt.xlabel('Частота, МГц') +plt.ylabel('Радианы') +plt.grid(True) + +# Действительная и мнимая части +plt.subplot(3, 2, 5) +plt.stem(freqs, real_part) +plt.title('Действительная часть $Re(p_T(f_n))$') +plt.xlabel('Частота, МГц') +plt.grid(True) + +plt.subplot(3, 2, 6) +plt.stem(freqs, imag_part) +plt.title('Мнимая часть $Im(p_T(f_n))$') +plt.xlabel('Частота, МГц') +plt.grid(True) + +plt.tight_layout() +for i, ax in enumerate(fig.axes): + bbox = ax.get_position() + print(f"График {i+1}: [{bbox.x0:.4f}, {bbox.y0:.4f}, {bbox.width:.4f}, {bbox.height:.4f}]") +plt.savefig('fig_1.png', dpi=300) diff --git a/code/3.m b/code/3.m new file mode 100644 index 0000000..1662de8 --- /dev/null +++ b/code/3.m @@ -0,0 +1,95 @@ +clear; clc; close all; + +blue_color = [0.12, 0.47, 0.71]; +red_color = [0.84, 0.15, 0.16]; +set(0, 'DefaultFigureColor', 'w'); +set(0, 'DefaultAxesColor', 'w'); +set(0, 'DefaultTextFontName', 'DejaVu Sans'); +set(0, 'DefaultAxesFontName', 'DejaVu Sans'); +set(0, 'DefaultAxesFontSize', 10); +set(0, 'DefaultTextFontSize', 12); +set(0, 'DefaultAxesLineWidth', 1.5); +set(0, 'DefaultLineLineWidth', 1.7); +set(0, 'DefaultTextFontWeight', 'normal'); +set(0, 'DefaultAxesTitleFontWeight', 'normal'); +set(0, 'DefaultAxesTickDir', 'in'); +set(0, 'DefaultAxesTickLength', [0.015 0.025]); + +set(0, 'DefaultAxesXGrid', 'on'); +set(0, 'DefaultAxesYGrid', 'on'); +set(0, 'DefaultAxesGridLineStyle', '-'); +set(0, 'DefaultAxesGridColor', [0.6 0.6 0.6]); +set(0, 'DefaultAxesGridAlpha', 0.5); + +a0 = 0.1; +f0 = 2.0; +w0 = 2 * pi * f0; +T = 0.5; +N = 16; +h = T / N; +fs = 1 / h; + +l_indices = 0:1:N; +t_l = l_indices * h; +p_l = 2 * a0 * sin(3*w0 * t_l) + a0 * cos(5 * w0 * t_l); +p_T_discrete = zeros(N); + +half_N = N / 2; +for n = 0 : half_N + sum_val = 0; + for l = 0:(N-1) + exponent = -1i * 2 * pi * n * l / N; + sum_val = sum_val + p_l(l+1) * exp(exponent); + end + p_T_discrete(n+1) = sum_val / N; +end + +for n = half_N+1:(N-1) + p_T_discrete(n) = conj(p_T_discrete((N - n) + 1)); +end + +n_extended = -N:1:N; +freqs_extended = n_extended .* (fs / N); +p_T_extended = p_T_discrete(mod(n_extended, N) + 1); +p_T_extended(abs(p_T_extended) < 1e-10) = 0; + +ampliude_ext = abs(p_T_extended); +phase_ext = angle(p_T_extended); +real_part_ext = real(p_T_extended); +imag_part_ext = imag(p_T_extended); + + +figure('Position', [110, 55, 1600, 1300]); + +axes('Position', [0.0628, 0.5508, 0.4298, 0.4129]); +stem(freqs_extended, ampliude_ext, 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +hold on; +y_bounds = ylim; +plot([fs/2, fs/2], y_bounds, '--r', 'LineWidth', 1.7, 'Color', red_color); +plot([-fs/2, -fs/2], y_bounds, '--r', 'LineWidth', 1.7, 'Color', red_color); +hold off; +title('\bf a) \rm Амплитудный спектр ДПФ'); +xlabel('Частота, МГц'); +ylabel('Амплитуда, МПа'); +grid on; + +axes('Position', [0.5577, 0.5508, 0.4298, 0.4129]); +stem(freqs_extended, phase_ext, 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color) +title('\bf b) \rm Фазовый спектр ДПФ'); +xlabel('Частота, МГц'); +ylabel('Фаза, рад'); +grid on; + +axes('Position', [0.0628, 0.0583, 0.4298, 0.4129]); +stem(freqs_extended, real_part_ext, 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color) +title('\bf c) \rm Действительная часть ДПФ'); +xlabel('Частота, МГц'); +grid on; + +axes('Position', [0.5577, 0.0583, 0.4298, 0.4129]); +stem(freqs_extended, imag_part_ext, 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color) +title('\bf d) \rm Мнимая часть ДПФ'); +xlabel('Частота, МГц'); +grid on; + +print(gcf, 'fig_3_m.png', '-dpng', '-r300'); diff --git a/code/3.py b/code/3.py new file mode 100644 index 0000000..9fefc92 --- /dev/null +++ b/code/3.py @@ -0,0 +1,79 @@ +import numpy as np +import matplotlib.pyplot as plt + + +a0 = 0.1 # MPa +f0 = 2.0 # MHz +w0 = 2 * np.pi * f0 +T = 0.5 # мкс +N = 16 # points +h = T / N # 0.03125 microseconds +fs = 1 / h # 32 MHz + + +l_indices = np.arange(N) +t_l = l_indices * h +p_l = 2 * a0 * np.sin(3 * w0 * t_l) + a0 * np.cos(5 * w0 * t_l) + +p_T_discrete = np.zeros(N, dtype=complex) + +half_N = N // 2 +for n in range(half_N + 1): + sum_val = 0j + for l in range(N): + exponent = -1j * 2 * np.pi * n * l / N + sum_val += p_l[l] * np.exp(exponent) + p_T_discrete[n] = sum_val / N + + +for n in range(half_N + 1, N): + p_T_discrete[n] = np.conj(p_T_discrete[N - n]) + +n_extended = np.arange(-N, N + 1) +freqs_extended = n_extended * (fs / N) +p_T_extended = p_T_discrete[n_extended % N] + +p_T_extended[np.abs(p_T_extended) < 1e-10] = 0 + +amplitude_ext = np.abs(p_T_extended) +phase_ext = np.angle(p_T_extended) +real_part_ext = np.real(p_T_extended) +imag_part_ext = np.imag(p_T_extended) + + +fig = plt.figure(figsize=(12, 10)) + +plt.subplot(2, 2, 1) +plt.stem(freqs_extended, amplitude_ext, basefmt=" ") +plt.title('Амплитудный спектр ДПФ') +plt.xlabel('Частота, МГц') +plt.ylabel('Амплитуда, МПа') +plt.axvline(x=fs/2, color='r', linestyle='--', alpha=0.5) +plt.axvline(x=-fs/2, color='r', linestyle='--', alpha=0.5) +plt.grid(True) + +plt.subplot(2, 2, 2) +plt.stem(freqs_extended, phase_ext, basefmt=" ") +plt.title('Фазовый спектр ДПФ') +plt.xlabel('Частота, МГц') +plt.ylabel('Фаза, рад') +plt.grid(True) + +plt.subplot(2, 2, 3) +plt.stem(freqs_extended, real_part_ext, basefmt=" ") +plt.title('Действительная часть ДПФ') +plt.xlabel('Частота, МГц') +plt.grid(True) + +plt.subplot(2, 2, 4) +plt.stem(freqs_extended, imag_part_ext, basefmt=" ") +plt.title('Мнимая часть ДПФ') +plt.xlabel('Частота, МГц') +plt.grid(True) + +plt.tight_layout() +print(fig.properties) +for i, ax in enumerate(fig.axes): + bbox = ax.get_position() + print(f"График {i+1}: [{bbox.x0:.4f}, {bbox.y0:.4f}, {bbox.width:.4f}, {bbox.height:.4f}]") +plt.savefig('fig_3.png', dpi=300) diff --git a/code/4.m b/code/4.m new file mode 100644 index 0000000..d6345d3 --- /dev/null +++ b/code/4.m @@ -0,0 +1,120 @@ +clear; clc; close all; + +blue_color = [0.12, 0.47, 0.71]; +red_color = [0.84, 0.15, 0.16]; +set(0, 'DefaultFigureColor', 'w'); +set(0, 'DefaultAxesColor', 'w'); +set(0, 'DefaultTextFontName', 'DejaVu Sans'); +set(0, 'DefaultAxesFontName', 'DejaVu Sans'); +set(0, 'DefaultAxesFontSize', 10); +set(0, 'DefaultTextFontSize', 12); +set(0, 'DefaultAxesLineWidth', 1.5); +set(0, 'DefaultLineLineWidth', 1.7); +set(0, 'DefaultTextFontWeight', 'normal'); +set(0, 'DefaultAxesTitleFontWeight', 'normal'); +set(0, 'DefaultAxesTickDir', 'in'); +set(0, 'DefaultAxesTickLength', [0.015 0.025]); + +set(0, 'DefaultAxesXGrid', 'on'); +set(0, 'DefaultAxesYGrid', 'on'); +set(0, 'DefaultAxesGridLineStyle', '-'); +set(0, 'DefaultAxesGridColor', [0.6 0.6 0.6]); +set(0, 'DefaultAxesGridAlpha', 0.5); + + +a0 = 0.1; +f0 = 2.0; +w0 = 2 * pi * f0; +T = 0.5; +N = 16; +h = T / N; + +l_indices = 0 : (N - 1); +t_l = l_indices * h; +p_l = 2 * a0 * sin(3 * w0 * t_l) + a0 * cos(5 * w0 * t_l); + +p_T_discrete = zeros(1, N); +for n = 0 : (N - 1) + sum_val = 0; + for l = 0 : (N - 1) + exponent = -1i * 2 * pi * n * l / N; + sum_val = sum_val + p_l(l + 1) * exp(exponent); + end + p_T_discrete(n + 1) = sum_val / N; +end + +p_restored = zeros(1, N); +for l = 0 : (N - 1) + sum_val = 0; + for n = 0 : (N - 1) + exponent = 1i * 2 * pi * n * l / N; + sum_val = sum_val + p_T_discrete(n + 1) * exp(exponent); + end + p_restored(l + 1) = sum_val; +end + +p_restored_real = real(p_restored); + +h_fine = h / 10; + +t_fine = 0 : h_fine : (T - h_fine); + +p_interpolated = zeros(1, length(t_fine)); + +for idx = 1 : length(t_fine) + t = t_fine(idx); + sum_sinc = 0; + for l = 0 : (N - 1) + sum_sinc = sum_sinc + p_l(l + 1) * sinc((t - l * h) / h); + end + p_interpolated(idx) = sum_sinc; +end + +p_analytical_fine = 2 * a0 * sin(3 * w0 * t_fine) + a0 * cos(5 * w0 * t_fine); + + +figure('Position', [110, 55, 1500, 1050]); + +axes('Position', [0.0781, 0.5651, 0.9069, 0.3890]); +scatter(t_l, p_l, 'filled', 150, 'o', 'LineWidth', 1.5, 'DisplayName', 'Исходный дискретный p(l)'); +hold on; +scatter(t_l, p_restored_real, 120, 'x', 'LineWidth', 1.5, 'DisplayName', 'Восстановленный ОДПФ'); +xlim([0, 0.45]); +hold off; + +title('\bf a) \rm Проверка обратимости прямого и обратного ДПФ'); +xlabel('Время t, мкс'); +ylabel('Амплитуда, МПа'); +grid on; +legend('Location', 'best'); + + +axes('Position', [0.0781, 0.0747, 0.9069, 0.3890]); +plot(t_l, p_l, 'ko', 'MarkerSize', 8, 'MarkerFaceColor', 'k', 'DisplayName', 'Узлы дискретизации p(l)'); +xlim([0, 0.45]); +hold on; +plot(t_fine, p_interpolated, 'r-', 'LineWidth', 2, 'DisplayName', 'Интерполяция (sinc)'); +plot(t_fine, p_analytical_fine, '--', 'Color', [0.5, 0.6, 1.0], 'LineWidth', 4, 'DisplayName', 'Аналитический сигнал'); +hold off; +title('\bf b) \rm Восстановление непрерывного сигнала (Теорема Котельникова)'); +xlabel('Время t, мкс'); +ylabel('Амплитуда, МПа'); +grid on; +legend('Location', 'best'); + + +print(gcf, 'fig_4_m.png', '-dpng', '-r300'); + +% evaluate absolute error +error_manual = p_restored_real - p_l; +max_err_manual = max(abs(error_manual)); +figure('Position', [150, 150, 900, 800]); +stem(t_l, error_manual, 'filled', 'Color', red_color, 'MarkerFaceColor', red_color, 'LineWidth', 1.5); +hold on; +plot(xlim, [0 0], 'k-', 'LineWidth', 1); +hold off; +title('Абсолютная ошибка восстановления точек для собственной реализации'); +xlabel('Время \itt\rm, мкс'); +ylabel('\Delta\itp\rm, МПа'); +grid on; +print(gcf, 'fig_4_error.png', '-dpng', '-r300'); diff --git a/code/4.py b/code/4.py new file mode 100644 index 0000000..ad9ee9a --- /dev/null +++ b/code/4.py @@ -0,0 +1,78 @@ +import numpy as np +import matplotlib.pyplot as plt + + +a0 = 0.1 +f0 = 2.0 +w0 = 2 * np.pi * f0 +T = 0.5 +N = 16 +h = T / N + +l_indices = np.arange(N) +t_l = l_indices * h +p_l = 2 * a0 * np.sin(3 * w0 * t_l) + a0 * np.cos(5 * w0 * t_l) + +p_T_discrete = np.zeros(N, dtype=complex) +for n in range(N): + sum_val = 0j + for l in range(N): + exponent = -1j * 2 * np.pi * n * l / N + sum_val += p_l[l] * np.exp(exponent) + p_T_discrete[n] = sum_val / N + + +p_restored = np.zeros(N, dtype=complex) +for l in range(N): + sum_val = 0j + for n in range(N): + exponent = 1j * 2 * np.pi * n * l / N + sum_val += p_T_discrete[n] * np.exp(exponent) + p_restored[l] = sum_val + + +p_restored_real = np.real(p_restored) + + +h_fine = h / 10 +t_fine = np.arange(0, T, h_fine) +p_interpolated = np.zeros(len(t_fine)) + +for idx, t in enumerate(t_fine): + sum_sinc = 0 + for l in range(N): + sum_sinc += p_l[l] * np.sinc((t - l * h) / h) + p_interpolated[idx] = sum_sinc + + +p_analytical_fine = 2 * a0 * np.sin(3 * w0 * t_fine) + a0 * np.cos(5 * w0 * t_fine) + + +fig = plt.figure(figsize=(10, 8)) + +plt.subplot(2, 1, 1) +plt.plot(t_l, p_l, 'bo-', markersize=8, label='Исходный дискретный $p(l)$') +plt.plot(t_l, p_restored_real, 'rX', markersize=6, label='Восстановленный ОДПФ') +plt.title('Проверка обратимости прямого и обратного ДПФ') +plt.xlabel('Время $t$, мкс') +plt.ylabel('Амплитуда, МПа') +plt.grid(True) +plt.legend() + +plt.subplot(2, 1, 2) +plt.plot(t_l, p_l, 'ko', markersize=6, label='Узлы дискретизации $p(l)$') +plt.plot(t_fine, p_interpolated, 'r-', linewidth=2, label='Интерполяция (sinc)') +plt.plot(t_fine, p_analytical_fine, 'b--', alpha=0.5, linewidth=4, label='Аналитический сигнал') + +plt.title('Восстановление непрерывного сигнала (Теорема Котельникова)') +plt.xlabel('Время $t$, мкс') +plt.ylabel('Амплитуда, МПа') +plt.grid(True) +plt.legend() + +plt.tight_layout() +print(fig.properties) +for i, ax in enumerate(fig.axes): + bbox = ax.get_position() + print(f"График {i+1}: [{bbox.x0:.4f}, {bbox.y0:.4f}, {bbox.width:.4f}, {bbox.height:.4f}]") +plt.savefig('fig_4.png', dpi=300) \ No newline at end of file diff --git a/code/5.m b/code/5.m new file mode 100644 index 0000000..2a33e58 --- /dev/null +++ b/code/5.m @@ -0,0 +1,83 @@ +clear; clc; close all; + +blue_color = [0.12, 0.47, 0.71]; +red_color = [0.84, 0.15, 0.16]; +set(0, 'DefaultFigureColor', 'w'); +set(0, 'DefaultAxesColor', 'w'); +set(0, 'DefaultTextFontName', 'DejaVu Sans'); +set(0, 'DefaultAxesFontName', 'DejaVu Sans'); +set(0, 'DefaultAxesFontSize', 10); +set(0, 'DefaultTextFontSize', 12); +set(0, 'DefaultAxesLineWidth', 1.5); +set(0, 'DefaultLineLineWidth', 1.7); +set(0, 'DefaultTextFontWeight', 'normal'); +set(0, 'DefaultAxesTitleFontWeight', 'normal'); +set(0, 'DefaultAxesTickDir', 'in'); +set(0, 'DefaultAxesTickLength', [0.015 0.025]); + +set(0, 'DefaultAxesXGrid', 'on'); +set(0, 'DefaultAxesYGrid', 'on'); +set(0, 'DefaultAxesGridLineStyle', '-'); +set(0, 'DefaultAxesGridColor', [0.6 0.6 0.6]); +set(0, 'DefaultAxesGridAlpha', 0.5); + + +a0 = 0.1; +f0 = 2.0; +w0 = 2 * pi * f0; +T = 0.5; +N = 16; +h = T / N; +l_indices = 0 : (N - 1); +t_l = l_indices * h; +p_l = 2 * a0 * sin(3 * w0 * t_l) + a0 * cos(5 * w0 * t_l); +P_fft_raw = fft(p_l); +P_fft_normalized = P_fft_raw / N; +amplitude_fft = abs(P_fft_normalized); +fs = 1 / h; +freqs_fft = (0 : (N - 1)) * (fs / N); +p_ifft = ifft(P_fft_raw); +p_ifft_real = real(p_ifft); + + +figure('Position', [110, 55, 1600, 1300]); +axes('Position', [0.0650, 0.5513, 0.9225, 0.4114]); +half_N = (N / 2) + 1; +stem(freqs_fft(1:half_N), amplitude_fft(1:half_N), 'filled', 'b', 'DisplayName', 'Спектр ДПФ'); +hold on; +plot(6.0, 0.1, 'ro', 'MarkerFaceColor', 'r', 'DisplayName', 'Аналитическая амплитуда (3f_0)'); +plot(10.0, 0.05, 'ro', 'MarkerFaceColor', 'r', 'DisplayName', 'Аналитическая амплитуда (5f_0)'); +hold off; +title('\bf a) \rm Амплитудный спектр, полученный с помощью БПФ (нормированный на N)'); +xlabel('Частота, МГц'); +ylabel('Амплитуда, МПа'); +grid on; +legend('Location', 'northeast'); + +axes('Position', [0.0650, 0.0598, 0.9225, 0.4114]); +scatter(t_l, p_l, 100, 'k-', 'LineWidth', 2, 'DisplayName', 'Исходный сигнал p(l)'); +hold on; +scatter(t_l, p_ifft_real, 120, 'x', 'LineWidth', 2, 'DisplayName', 'Восстановленный ifft'); +hold off; +title('\bf b) \rm Проверка обратимости прямого и обратного БПФ (FFT / IFFT)'); +xlabel('Время \itt\rm, мкс'); +ylabel('Амплитуда, МПа'); +grid on; +legend('Location', 'northeast'); +print(gcf, 'fig_5_m.png', '-dpng', '-r300'); + +% evaluate absolute error +error_fft = p_ifft_real - p_l; +max_err_fft = max(abs(error_fft)); +figure('Position', [150, 150, 900, 800]); +stem(t_l, error_fft, 'filled', 'Color', red_color, 'MarkerFaceColor', red_color, 'LineWidth', 1.5); +hold on; +plot(xlim, [0 0], 'k-', 'LineWidth', 1); +hold off; +title('Абсолютная ошибка восстановления точек для встроенной реализации'); +xlabel('Время \itt\rm, мкс'); +ylabel('\Delta\itp\rm, МПа'); +grid on; +text_str2 = sprintf('Max Error = %.2e', max_err_fft); +legend(text_str2, 'Location', 'best'); +print(gcf, 'fig_5_error.png', '-dpng', '-r300'); diff --git a/code/5.py b/code/5.py new file mode 100644 index 0000000..b3ce702 --- /dev/null +++ b/code/5.py @@ -0,0 +1,62 @@ +import numpy as np +import matplotlib.pyplot as plt + + +a0 = 0.1 +f0 = 2.0 +w0 = 2 * np.pi * f0 +T = 0.5 +N = 16 +h = T / N + +l_indices = np.arange(N) +t_l = l_indices * h +p_l = 2 * a0 * np.sin(3 * w0 * t_l) + a0 * np.cos(5 * w0 * t_l) + + +P_fft_raw = np.fft.fft(p_l) + + +P_fft_normalized = P_fft_raw / N +amplitude_fft = np.abs(P_fft_normalized) + + +fs = 1 / h +freqs_fft = np.arange(N) * (fs / N) + + +p_ifft = np.fft.ifft(P_fft_raw) +p_ifft_real = np.real(p_ifft) + + +fig = plt.figure(figsize=(12, 10)) + +plt.subplot(2, 1, 1) + +half_N = N // 2 + 1 +plt.stem(freqs_fft[:half_N], amplitude_fft[:half_N], basefmt=" ",) + +plt.plot(6.0, 0.1, 'ro', label='Аналитическая амплитуда (3f0)') +plt.plot(10.0, 0.05, 'ro', label='Аналитическая амплитуда (5f0)') + +plt.title('Амплитудный спектр, полученный с помощью БПФ (нормированный на $N$)') +plt.xlabel('Частота, МГц') +plt.ylabel('Амплитуда, МПа') +plt.grid(True) +plt.legend() + +plt.subplot(2, 1, 2) +plt.plot(t_l, p_l, 'k-', linewidth=2, label='Исходный сигнал $p(l)$') +plt.plot(t_l, p_ifft_real, 'r--', linewidth=2, label='Восстановленный $ifft$') +plt.title('Проверка обратимости прямого и обратного БПФ (FFT / IFFT)') +plt.xlabel('Время $t$, мкс') +plt.ylabel('Амплитуда, МПа') +plt.grid(True) +plt.legend() + +plt.tight_layout() +print(fig.properties) +for i, ax in enumerate(fig.axes): + bbox = ax.get_position() + print(f"График {i+1}: [{bbox.x0:.4f}, {bbox.y0:.4f}, {bbox.width:.4f}, {bbox.height:.4f}]") +plt.savefig('fig_5.png', dpi=300) diff --git a/code/6.m b/code/6.m new file mode 100644 index 0000000..71ac3ae --- /dev/null +++ b/code/6.m @@ -0,0 +1,160 @@ +clear; clc; close all; + +blue_color = [0.12, 0.47, 0.71]; +red_color = [0.84, 0.15, 0.16]; +set(0, 'DefaultFigureColor', 'w'); +set(0, 'DefaultAxesColor', 'w'); +set(0, 'DefaultTextFontName', 'DejaVu Sans'); +set(0, 'DefaultAxesFontName', 'DejaVu Sans'); +set(0, 'DefaultAxesFontSize', 10); +set(0, 'DefaultTextFontSize', 12); +set(0, 'DefaultAxesLineWidth', 1.5); +set(0, 'DefaultLineLineWidth', 1.7); +set(0, 'DefaultTextFontWeight', 'normal'); +set(0, 'DefaultAxesTitleFontWeight', 'normal'); +set(0, 'DefaultAxesTickDir', 'in'); +set(0, 'DefaultAxesTickLength', [0.015 0.025]); + +set(0, 'DefaultAxesXGrid', 'on'); +set(0, 'DefaultAxesYGrid', 'on'); +set(0, 'DefaultAxesGridLineStyle', '-'); +set(0, 'DefaultAxesGridColor', [0.6 0.6 0.6]); +set(0, 'DefaultAxesGridAlpha', 0.5); + + +a0 = 0.1; +f0 = 2.0; +w0 = 2 * pi * f0; +T = 0.5; +N = 16; +h = T / N; +fs = 1 / h; +l_indices = 0 : (N - 1); +t_l = l_indices * h; +p_l = 2 * a0 * sin(3 * w0 * t_l) + a0 * cos(5 * w0 * t_l); +freqs_f0_multipliers = -6 : 6; +freqs = freqs_f0_multipliers * f0; +coeffs_deriv = zeros(1, length(freqs)); +coeffs_deriv(freqs_f0_multipliers == 3) = 0.6 * pi * f0; +coeffs_deriv(freqs_f0_multipliers == -3) = 0.6 * pi * f0; +coeffs_deriv(freqs_f0_multipliers == 5) = 1i * 0.5 * pi * f0; +coeffs_deriv(freqs_f0_multipliers == -5) = -1i * 0.5 * pi * f0; + + +figure('Position', [110, 55, 1600, 1300]); + +axes('Position', [0.0391, 0.5390, 0.4498, 0.4156]); +stem(freqs, abs(coeffs_deriv), 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +title('\bf a) \rm Амплитудный спектр dp/dt'); + +axes('Position', [0.5377, 0.5390, 0.4498, 0.4156]); +stem(freqs, angle(coeffs_deriv), 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +title('\bf b) \rm Фазовый спектр dp/dt'); + +axes('Position', [0.0391, 0.0484, 0.4498, 0.4156]); +stem(freqs, real(coeffs_deriv), 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +title('\bf c) \rm Действительная часть dp/dt'); + +axes('Position', [0.5377, 0.0484, 0.4498, 0.4156]); +stem(freqs, imag(coeffs_deriv), 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +title('\bf d) \rm Мнимая часть dp/dt'); + +print(gcf, 'fig_add_1_m.png', '-dpng', '-r300'); + +p_T_discrete = zeros(1, N); +for n = 0 : (N - 1) + sum_val = 0; + for l = 0 : (N - 1) + exponent = -1i * 2 * pi * n * l / N; + sum_val = sum_val + p_l(l + 1) * exp(exponent); + end + p_T_discrete(n + 1) = sum_val / N; +end + +freqs_mapped = zeros(1, N); +for n = 0 : (N - 1) + if n <= N/2 + freqs_mapped(n + 1) = n * (fs / N); + else + freqs_mapped(n + 1) = (n - N) * (fs / N); + end +end + +omega_mapped = 2 * pi * freqs_mapped; + +p_T_deriv_discrete = p_T_discrete .* (1i * omega_mapped); + +n_extended = -N : N; +freqs_extended = n_extended * (fs / N); +p_T_deriv_extended = p_T_deriv_discrete(mod(n_extended, N) + 1); +p_T_deriv_extended(abs(p_T_deriv_extended) < 1e-10) = 0; + +figure('Position', [110, 55, 1500, 1050]); + +axes('Position', [0.0391, 0.5390, 0.4498, 0.4141]); +stem(freqs_extended, abs(p_T_deriv_extended), 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +hold on; +y_bounds = ylim; +plot([fs/2, fs/2], y_bounds, '--', 'Color', red_color, 'LineWidth', 1.5); +plot([-fs/2, -fs/2], y_bounds, '--', 'Color', red_color, 'LineWidth', 1.5); +hold off; +title('\bf a) \rm Амплитуда P_discrete в интервале [-1/h, 1/h]'); + + +axes('Position', [0.5377, 0.5390, 0.4498, 0.4141]); +stem(freqs_extended, angle(p_T_deriv_extended), 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +title('\bf b) \rm Фаза P_{discrete}'); + + +axes('Position', [0.0391, 0.0484, 0.4498, 0.4141]); +stem(freqs_extended, real(p_T_deriv_extended), 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +title('\bf c) \rm Re(P_{discrete})'); + + +axes('Position', [0.5377, 0.0484, 0.4498, 0.4141]); +stem(freqs_extended, imag(p_T_deriv_extended), 'filled', 'Color', blue_color, 'MarkerFaceColor', blue_color); +title('\bf d) \rm Im(P_{discrete})'); + +print(gcf, 'fig_add_2_m.png', '-dpng', '-r300'); + +P_fft = fft(p_l); + +% there is no fftfreq in matlab, so these freqs came out from python script +freqs_fft = [0 : (N/2 - 1), -N/2 : -1] * (fs / N); +omega_fft = 2 * pi * freqs_fft; +P_deriv_fft = P_fft .* (1i * omega_fft); +p_deriv_ifft = ifft(P_deriv_fft); +p_deriv_ifft_real = real(p_deriv_ifft); +p_deriv_analytical = 0.6 * w0 * cos(3 * w0 * t_l) - 0.5 * w0 * sin(5 * w0 * t_l); + +figure('Position', [110, 55, 1500, 750]); +scatter(t_l, p_deriv_analytical, 100, 'k-', 'LineWidth', 3, 'DisplayName', 'Аналитическая dp/dt'); +hold on; +scatter(t_l, p_deriv_ifft_real, 120, 'x-', 'LineWidth', 2, 'DisplayName', 'Восстановленная ifft'); +hold off; + +title('Восстановление производной через обратное БПФ'); +xlabel('Время t, мкс'); +ylabel('Амплитуда, МПа/мкс'); +legend('Location', 'northeast'); + +print(gcf, 'fig_add_3_m.png', '-dpng', '-r300'); + + +% evaluate absolute error for derivative +error_deriv = p_deriv_ifft_real - p_deriv_analytical; +max_err_deriv = max(abs(error_deriv)); +figure('Position', [110, 55, 1500, 600]); +stem(t_l, error_deriv, 'filled', 'Color', red_color, 'MarkerFaceColor', red_color); + +hold on; +y_bounds = ylim; +plot(xlim, [0 0], 'k-', 'LineWidth', 1); +hold off; + +title('Абсолютная ошибка вычисления производной через БПФ'); +xlabel('Время t, мкс'); +ylabel('\Delta dp/dt, МПа/мкс'); +text_str = sprintf('Max Error = %.2e', max_err_deriv); +legend(text_str, 'Location', 'northeast'); +print(gcf, 'fig_add_6_error.png', '-dpng', '-r300'); diff --git a/code/6.py b/code/6.py new file mode 100644 index 0000000..f107e10 --- /dev/null +++ b/code/6.py @@ -0,0 +1,143 @@ +import numpy as np +import matplotlib.pyplot as plt + +a0 = 0.1 +f0 = 2.0 +w0 = 2 * np.pi * f0 # 4*/pi +T = 0.5 +N = 16 +h = T / N +fs = 1 / h + + +l_indices = np.arange(N) +t_l = l_indices * h +p_l = 2 * a0 * np.sin(3 * w0 * t_l) + a0 * np.cos(5 * w0 * t_l) + +# 1 +freqs_f0_multipliers = np.arange(-6, 7) +freqs = freqs_f0_multipliers * f0 + + +coeffs_deriv = np.zeros(len(freqs), dtype=complex) +coeffs_deriv[freqs_f0_multipliers == 3] = 0.6 * np.pi * f0 # 1.2 * /pi +coeffs_deriv[freqs_f0_multipliers == -3] = 0.6 * np.pi * f0 +coeffs_deriv[freqs_f0_multipliers == 5] = 1j * 0.5 * np.pi * f0 # 1j * /pi +coeffs_deriv[freqs_f0_multipliers == -5] = -1j * 0.5 * np.pi * f0 + +fig = plt.figure(figsize=(12, 8)) +plt.subplot(2, 2, 1) +plt.stem(freqs, np.abs(coeffs_deriv), basefmt=" ") +plt.title('Амплитудный спектр dp/dt') +plt.grid(True) + +plt.subplot(2, 2, 2) +plt.stem(freqs, np.angle(coeffs_deriv), basefmt=" ") +plt.title('Фазовый спектр dp/dt') +plt.grid(True) + +plt.subplot(2, 2, 3) +plt.stem(freqs, np.real(coeffs_deriv), basefmt=" ") +plt.title('Действительная часть dp/dt') +plt.grid(True) + +plt.subplot(2, 2, 4) +plt.stem(freqs, np.imag(coeffs_deriv), basefmt=" ") +plt.title('Мнимая часть dp/dt') +plt.grid(True) +plt.tight_layout() +print(fig.properties) +for i, ax in enumerate(fig.axes): + bbox = ax.get_position() + print(f"График {i+1}: [{bbox.x0:.4f}, {bbox.y0:.4f}, {bbox.width:.4f}, {bbox.height:.4f}]") +plt.savefig('fig_add_1.png', dpi=300) +plt.close() + + +# 2 +p_T_discrete = np.zeros(N, dtype=complex) +for n in range(N): + sum_val = 0j + for l in range(N): + exponent = -1j * 2 * np.pi * n * l / N + sum_val += p_l[l] * np.exp(exponent) + p_T_discrete[n] = sum_val / N + + +freqs_mapped = np.zeros(N) +for n in range(N): + if n <= N//2: + freqs_mapped[n] = n * (fs / N) + else: + freqs_mapped[n] = (n - N) * (fs / N) + +omega_mapped = 2 * np.pi * freqs_mapped + +p_T_deriv_discrete = p_T_discrete * (1j * omega_mapped) + +n_extended = np.arange(-N, N + 1) +freqs_extended = n_extended * (fs / N) +p_T_deriv_extended = p_T_deriv_discrete[n_extended % N] +p_T_deriv_extended[np.abs(p_T_deriv_extended) < 1e-10] = 0 # get rid of noise + +fig = plt.figure(figsize=(12, 8)) +plt.subplot(2, 2, 1) +plt.stem(freqs_extended, np.abs(p_T_deriv_extended), basefmt=" ") +plt.title('Амплитуда $P\'_{discrete}$ в интервале $[-1/h, 1/h]$') +plt.axvline(x=fs/2, color='r', linestyle='--', alpha=0.5) +plt.axvline(x=-fs/2, color='r', linestyle='--', alpha=0.5) +plt.grid(True) + +plt.subplot(2, 2, 2) +plt.stem(freqs_extended, np.angle(p_T_deriv_extended), basefmt=" ") +plt.title('Фаза $P\'_{discrete}$') +plt.grid(True) + +plt.subplot(2, 2, 3) +plt.stem(freqs_extended, np.real(p_T_deriv_extended), basefmt=" ") +plt.title('Re($P\'_{discrete}$)') +plt.grid(True) + +plt.subplot(2, 2, 4) +plt.stem(freqs_extended, np.imag(p_T_deriv_extended), basefmt=" ") +plt.title('Im($P\'_{discrete}$)') +plt.grid(True) +plt.tight_layout() +print(fig.properties) +for i, ax in enumerate(fig.axes): + bbox = ax.get_position() + print(f"График {i+1}: [{bbox.x0:.4f}, {bbox.y0:.4f}, {bbox.width:.4f}, {bbox.height:.4f}]") +plt.savefig('fig_add_2.png', dpi=300) +plt.close() + + +# 3* +P_fft = np.fft.fft(p_l) + +freqs_fft = np.fft.fftfreq(N, d=h) +print(freqs_fft) +omega_fft = 2 * np.pi * freqs_fft + +P_deriv_fft = P_fft * (1j * omega_fft) + +p_deriv_ifft = np.fft.ifft(P_deriv_fft) +p_deriv_ifft_real = np.real(p_deriv_ifft) + +# dp/dt = 0.6 * w0 * cos(3*w0*t) - 0.5 * w0 * sin(5*w0*t) +p_deriv_analytical = 0.6 * w0 * np.cos(3 * w0 * t_l) - 0.5 * w0 * np.sin(5 * w0 * t_l) + +fig = plt.figure(figsize=(10, 5)) +plt.plot(t_l, p_deriv_analytical, 'k-', linewidth=3, label='Аналитическая $dp/dt$') +plt.plot(t_l, p_deriv_ifft_real, 'r--', marker='o', markersize=6, linewidth=2, label='Восстановленная ifft( P * iw )') +plt.title('Восстановление производной через обратное БПФ') +plt.xlabel('Время $t$, мкс') +plt.ylabel('Амплитуда, МПа/мкс') +plt.legend() +plt.grid(True) +plt.tight_layout() +print(fig.properties) +for i, ax in enumerate(fig.axes): + bbox = ax.get_position() + print(f"График {i+1}: [{bbox.x0:.4f}, {bbox.y0:.4f}, {bbox.width:.4f}, {bbox.height:.4f}]") +plt.savefig('fig_add_3.png', dpi=300) +plt.close() diff --git a/code/fig_1_m.png b/code/fig_1_m.png new file mode 100644 index 0000000..869a99e Binary files /dev/null and b/code/fig_1_m.png differ diff --git a/code/fig_3_m.png b/code/fig_3_m.png new file mode 100644 index 0000000..1ec306e Binary files /dev/null and b/code/fig_3_m.png differ diff --git a/code/fig_4_error.png b/code/fig_4_error.png new file mode 100644 index 0000000..903a532 Binary files /dev/null and b/code/fig_4_error.png differ diff --git a/code/fig_4_m.png b/code/fig_4_m.png new file mode 100644 index 0000000..d0a84fb Binary files /dev/null and b/code/fig_4_m.png differ diff --git a/code/fig_5_error.png b/code/fig_5_error.png new file mode 100644 index 0000000..129c78f Binary files /dev/null and b/code/fig_5_error.png differ diff --git a/code/fig_5_m.png b/code/fig_5_m.png new file mode 100644 index 0000000..fbb8ac5 Binary files /dev/null and b/code/fig_5_m.png differ diff --git a/code/fig_add_1_m.png b/code/fig_add_1_m.png new file mode 100644 index 0000000..8810983 Binary files /dev/null and b/code/fig_add_1_m.png differ diff --git a/code/fig_add_2_m.png b/code/fig_add_2_m.png new file mode 100644 index 0000000..56d07a3 Binary files /dev/null and b/code/fig_add_2_m.png differ diff --git a/code/fig_add_3_m.png b/code/fig_add_3_m.png new file mode 100644 index 0000000..8de9fa2 Binary files /dev/null and b/code/fig_add_3_m.png differ diff --git a/code/fig_add_6_error.png b/code/fig_add_6_error.png new file mode 100644 index 0000000..cfed51c Binary files /dev/null and b/code/fig_add_6_error.png differ diff --git a/code/octave-workspace b/code/octave-workspace new file mode 100644 index 0000000..16479f0 Binary files /dev/null and b/code/octave-workspace differ diff --git a/fig_3.png b/fig_3.png new file mode 100644 index 0000000..359ca54 Binary files /dev/null and b/fig_3.png differ diff --git a/fig_4.png b/fig_4.png new file mode 100644 index 0000000..03511cd Binary files /dev/null and b/fig_4.png differ diff --git a/fig_5.png b/fig_5.png new file mode 100644 index 0000000..53c7270 Binary files /dev/null and b/fig_5.png differ diff --git a/fig_add_1.png b/fig_add_1.png new file mode 100644 index 0000000..e8333c6 Binary files /dev/null and b/fig_add_1.png differ diff --git a/fig_add_2.png b/fig_add_2.png new file mode 100644 index 0000000..68bcaeb Binary files /dev/null and b/fig_add_2.png differ diff --git a/fig_add_3.png b/fig_add_3.png new file mode 100644 index 0000000..e53479a Binary files /dev/null and b/fig_add_3.png differ diff --git a/pic/fig_1.png b/pic/fig_1.png new file mode 100644 index 0000000..ef96b12 Binary files /dev/null and b/pic/fig_1.png differ diff --git a/pic/fig_1_m.png b/pic/fig_1_m.png new file mode 100644 index 0000000..39c88f5 Binary files /dev/null and b/pic/fig_1_m.png differ diff --git a/pic/fig_3.png b/pic/fig_3.png new file mode 100644 index 0000000..359ca54 Binary files /dev/null and b/pic/fig_3.png differ diff --git a/pic/fig_3_m.png b/pic/fig_3_m.png new file mode 100644 index 0000000..1ec306e Binary files /dev/null and b/pic/fig_3_m.png differ diff --git a/pic/fig_4.png b/pic/fig_4.png new file mode 100644 index 0000000..03511cd Binary files /dev/null and b/pic/fig_4.png differ diff --git a/pic/fig_4_error.png b/pic/fig_4_error.png new file mode 100644 index 0000000..903a532 Binary files /dev/null and b/pic/fig_4_error.png differ diff --git a/pic/fig_4_m.png b/pic/fig_4_m.png new file mode 100644 index 0000000..d0a84fb Binary files /dev/null and b/pic/fig_4_m.png differ diff --git a/pic/fig_5.png b/pic/fig_5.png new file mode 100644 index 0000000..53c7270 Binary files /dev/null and b/pic/fig_5.png differ diff --git a/pic/fig_5_error.png b/pic/fig_5_error.png new file mode 100644 index 0000000..129c78f Binary files /dev/null and b/pic/fig_5_error.png differ diff --git a/pic/fig_5_m.png b/pic/fig_5_m.png new file mode 100644 index 0000000..fbb8ac5 Binary files /dev/null and b/pic/fig_5_m.png differ diff --git a/pic/fig_add_1.png b/pic/fig_add_1.png new file mode 100644 index 0000000..e8333c6 Binary files /dev/null and b/pic/fig_add_1.png differ diff --git a/pic/fig_add_1_m.png b/pic/fig_add_1_m.png new file mode 100644 index 0000000..8810983 Binary files /dev/null and b/pic/fig_add_1_m.png differ diff --git a/pic/fig_add_2.png b/pic/fig_add_2.png new file mode 100644 index 0000000..68bcaeb Binary files /dev/null and b/pic/fig_add_2.png differ diff --git a/pic/fig_add_2_m.png b/pic/fig_add_2_m.png new file mode 100644 index 0000000..56d07a3 Binary files /dev/null and b/pic/fig_add_2_m.png differ diff --git a/pic/fig_add_3_m.png b/pic/fig_add_3_m.png new file mode 100644 index 0000000..8de9fa2 Binary files /dev/null and b/pic/fig_add_3_m.png differ diff --git a/pic/fig_add_6_error.png b/pic/fig_add_6_error.png new file mode 100644 index 0000000..cfed51c Binary files /dev/null and b/pic/fig_add_6_error.png differ diff --git a/slides/.DS_Store b/slides/.DS_Store new file mode 100644 index 0000000..431ac28 Binary files /dev/null and b/slides/.DS_Store differ diff --git a/slides/presentation.md b/slides/presentation.md new file mode 100644 index 0000000..9a1676e --- /dev/null +++ b/slides/presentation.md @@ -0,0 +1,303 @@ +--- +marp: false +theme: default +paginate: true +backgroundColor: #ffffff +math: mathjax +style: | + section { + font-size: 26px; /* Чуть уменьшили шрифт */ + justify-content: start; + } + h1, h2 { + color: #1f77b4; + } + .lead { + display: flex; + flex-direction: column; + justify-content: center; + align-items: center; + text-align: center; + } + /* Делаем формулы чуть компактнее */ + p { + margin-bottom: 0.5em; + } +--- + + + + +# Дискретное преобразование Фурье + +**Самостоятельная работа. Вариант 1** +*Парфенов Михаил, 424 группа* + +--- + + +## 1. Постановка задачи + +Требуется исследовать акустический сигнал, описываемый бигармонической функцией: + +$$ +p(t) = 2a_0 \sin(3\omega_0 t) + a_0 \cos(5\omega_0 t) +$$ + +**Заданные параметры системы:** +* Амплитуда: $a_0 = 0.1$ МПа +* Базовая частота: $f_0 = 2$ МГц +* Круговая частота: $\omega_0 = 2\pi f_0$ + +**Цель работы:** +Реализация алгоритмов ДПФ/ОДПФ, интерполяция сигнала и спектральное дифференцирование. + +--- + +![bg right:65% 85%](../pic/fig_1_m.png) + + +## 2. Аналитический спектр + +**Гармоники:** $f_1 = 6$ МГц, $f_2 = 10$ МГц. +**Период:** $T = 1 / НОД = \mathbf{0.5 \text{ мкс}}$. + +**Амплитуды:** +* $\pm 6$ МГц: $\mp 0.1i$ +* $\pm 10$ МГц: $0.05$ + +**Симметрия:** +* $|P|$ и $\text{Re}$ — **четные**. +* $\arg$ и $\text{Im}$ — **нечетные**. + +--- + + + +## 3. Выбор параметров дискретизации + +Согласно **теореме Котельникова-Шеннона**: +$$ f_s > 2f_{max} \implies f_s > 20 \text{ МГц} $$ + +**Расчет сетки:** +* Шаг дискретизации: $h \le 0.05$ мкс. +* Окно дискретизации: $T = 0.5$ мкс. +* Минимальное число точек: $N^* = T/h = 10$. + +Для оптимизации алгоритма БПФ выбрано ближайшее число, являющееся степенью двойки: **$N = 16$**. +*(Фактический шаг: $h = 0.03125$ мкс, $f_s = 32$ МГц).* + +--- + +![bg right:65% 85%](../pic/fig_3_m.png) + + +## 4. Собственная реализация ДПФ + +Вычисление спектра: +$$ \overline{p_T(n)} = \frac{1}{N} \sum_{l=0}^{N-1} p(l) \exp\left(-i \frac{2\pi nl}{N} \right) $$ + +**Оптимизация:** +Рассчитывалась только первая половина спектра ($0 \dots N/2$). Вторая получена сопряжением: +$$ p_T(N - n) = p_T^*(n) $$ + +*Амплитуды в 2 раза меньше аналитических из-за двустороннего спектра.* + +--- + + +![bg right:65% 85%](../pic/fig_4_m.png) + + +## 5. Обратное ДПФ и Интерполяция + +**Восстановление:** +$$ p(l) = \sum_{n=0}^{N-1} \overline{p_T(n)} \exp\left( i\frac{2\pi nl}{N} \right) $$ + +**Теорема Котельникова:** +Реконструкция непрерывного сигнала: +$$ p(t) = \sum_{l=0}^{N-1} p(l)\operatorname{sinc}\left[ \frac{\pi}{h}(t - lh) \right] $$ + +*Ошибка восстановления $\sim 10^{-16}$.* + +--- + + +![bg right:65% 85%](../pic/fig_5_m.png) + + +## 6. Встроенные алгоритмы FFT + +Использованы функции `fft` и `ifft`. + +**Масштабирование:** +Встроенные алгоритмы вычисляют ненормированную сумму. Для получения физических амплитуд применен коэффициент $1/N$: +$$ P_{physical} = \frac{fft(p)}{N} $$ + +*Алгоритмы составляют обратимую пару (ошибка $\sim 10^{-17}$).* + +--- + + +## 7. Дополнительная задача: Дифференцирование + +Дифференцирование сигнала во временной области эквивалентно умножению его спектра на $i\omega$ в частотной: +$$ \frac{dp(t)}{dt} \Longleftrightarrow i\omega_n \cdot P(\omega_n) $$ + +**Алгоритм реализации:** +1. Прямое БПФ исходного сигнала. +2. Формирование вектора частот $\omega_n$ с учетом отрицательной полуоси. +3. Умножение: $P'_{fft} = P_{fft} \cdot (i\omega_n)$. +4. Обратное БПФ: $p' = real(ifft(P'_{fft}))$. + +--- + + +## Спектр производной + +Аналитически можно получить: +* Для $\pm 3f_0$: $C_{\pm 3} = \frac{0.6\omega_0}{2} = \mathbf{0.6\pi f_0}$ +* Для $\pm 5f_0$: $C_{\pm 5} = \mp i \frac{0.5\omega_0}{2} = \mathbf{\mp i 0.5\pi f_0}$ + +Что совпадает с численными результатами: + +![bg right:65% 85%](../pic/fig_add_1_m.png) + + +--- + +![bg right:60% contain](../pic/fig_add_3_m.png) + + + +## Результат дифференцирования + +Восстановленная производная $dp/dt$ совпадает с аналитически вычисленной производной. + +Максимальная ошибка восстановления составляет $\sim 10^{-14}$ МПа/мкс. + +--- + + +## Анализ точности восстановления + +$$ \Delta p(t) = p_{restored}(t) - p_{original}(t) $$ + +
+ + +
+ +--- + + + + +## Эмуляция стиля Matplotlib в MATLAB + +Дефолтные графики MATLAB (Octave) выглядят устаревшими. Чтобы отчет выглядел современно, стиль **Matplotlib** был воссоздан "с нуля" через глобальные настройки (аналог `plt.rcParams`). + +
+
+ +**Главные хаки:** +1. **Цветовая палитра:** Вручную заданы RGB-коды дефолтных питоновских цветов (Tableau 10). Знакомый синий `C0` = `[0.12, 0.47, 0.71]`. +2. **Шрифты:** Жирность заголовков принудительно отключена (`normal`), установлен питоновский шрифт `DejaVu Sans`. +3. **Стемы:** Магический параметр `'filled'` в функции `stem()` закрашивает пустые кружки. +4. **Базовая линия:** В Octave красная ось $Y=0$ реализуется отрисовкой дополнительной линии: `plot(xlim, [0 0], 'r')`. + +
+
+ +**Фрагмент конфигурации (MATLAB):** +```matlab +% Задаем цвета Python (C0 и C3) +mpl_blue = [0.12, 0.47, 0.71]; +mpl_red = [0.84, 0.15, 0.16]; + +% Настройки шрифтов и фона +set(0, 'DefaultFigureColor', 'w'); +set(0, 'DefaultTextFontName', 'DejaVu Sans'); +set(0, 'DefaultAxesTitleFontWeight', 'normal'); + +% Делаем сетку сплошной и полупрозрачной +set(0, 'DefaultAxesGridLineStyle', '-'); +set(0, 'DefaultAxesGridColor', [0.6 0.6 0.6]); + +% Рисуем красивый stem +stem(freqs, amp, 'filled', 'Color', mpl_blue); +``` +
+
+ +--- + + + +### Вопросы? + +*Примеры вопросов*: +1) когда отменят карантин на фф? +2) как сделать такую же классную презентацию? + +
+ +
+ +--- \ No newline at end of file diff --git a/slides/presentation.pdf b/slides/presentation.pdf new file mode 100644 index 0000000..dbf0532 Binary files /dev/null and b/slides/presentation.pdf differ diff --git a/slides/qr.png b/slides/qr.png new file mode 100644 index 0000000..1bee4ef Binary files /dev/null and b/slides/qr.png differ