PySDR学习之旅——2.频域

最新推荐文章于 2024-08-30 07:41:08 发布

彪嘞个彪

最新推荐文章于 2024-08-30 07:41:08 发布

阅读量255

点赞数

分类专栏： PySDR学习之旅文章标签：学习

原文链接：https://pysdr.org/content/frequency_domain.html

版权

PySDR学习之旅专栏收录该内容

1 篇文章 0 订阅

订阅专栏

本章使用 Python 示例介绍频域，并涵盖傅里叶级数、傅里叶变换、傅里叶性质、FFT、加窗和频谱图。

学习 DSP 和无线通信最酷的副作用之一是您还将学会在频域中思考。大多数人在频域工作的经验仅限于调整汽车音响系统上的低音/中音/高音旋钮。大多数人在频域中查看某些内容的经验仅限于查看音频均衡器，例如此剪辑：

在本章结束时，您将了解频域的真正含义，如何在时间和频率之间进行转换（以及这样做时会发生什么），以及我们将在 DSP 和 SDR 研究中使用的一些有趣原理。学完本教科书后，保证您将成为频域工作的大师！

首先，为什么我们喜欢在频域中查看信号？这里有两个示例信号，在时域和频域中显示。

正如您所看到的，在时域中它们看起来都像噪声，但在频域中我们可以看到不同的特征。一切事物都以其自然的形式存在于时域中；当我们对信号进行采样时，我们将在时域中对它们进行采样，因为您不能直接在频域中对信号进行采样。但有趣的事情通常发生在频域中。

傅立叶级数

频域的基础知识首先要理解任何信号都可以用正弦波相加来表示。当我们将信号分解为其复合正弦波时，我们将其称为傅里叶级数。以下是仅由两个正弦波组成的信号示例：

这是另一个例子；下面的红色曲线通过将 10 个正弦波相加来近似锯齿波。我们可以看到，这不是一个完美的重建——由于急剧的转变，需要无限数量的正弦波来重现这个锯齿波：

有些信号比其他信号需要更多的正弦波，有些信号需要无限的数量，尽管它们总是可以用有限的数量来近似。这是信号被分解为一系列正弦波的另一个示例：

要了解如何将信号分解为正弦波或正弦波，我们需要首先回顾一下正弦波的三个属性：

振幅
频率
阶段

振幅表示波的“强度”，而频率是每秒波的数量。相位用于表示正弦波如何随时间移动，范围从 0 到 360 度（或 0 到），但它必须相对于某事物才具有任何意义，例如具有相同频率的两个信号相差 30 度彼此的相位。

此时您可能已经意识到“信号”本质上只是一个函数，通常表示“随时间变化”（即 x 轴）。另一个容易记住的属性是周期，它是频率的倒数。正弦曲线的周期是波完成一个周期的时间量（以秒为单位）。因此，频率的单位是1/秒，即Hz。

当我们将信号分解为正弦波的总和时，每个正弦波都有一定的幅度、相位和频率。每个正弦波的幅度将告诉我们原始信号中存在的频率有多强。现在不要太担心相位，除了意识到 sin() 和 cos() 之间的唯一区别是相移（时间偏移）。

理解基本概念比求解傅立叶级数的实际方程更重要，但对于那些对方程感兴趣的人，我建议您参考 Wolfram 的简洁解释： https: //mathworld.wolfram.com/FourierSeries.html。

时频对

我们已经确定信号可以表示为正弦波，它具有多个属性。现在，让我们学习在频域中绘制信号。时域展示了信号如何随时间变化，而频域则展示了信号有多少位于哪些频率。X 轴不再是时间，而是频率。我们可以在时间和频率上绘制给定信号。让我们先看一些简单的例子。

以下是频率为 f 的正弦波在时域和频域中的样子：

时域看起来应该非常熟悉。这是一个振荡函数。不要担心它在周期的哪个点开始或持续多长时间。结论是信号具有单一频率，这就是为什么我们在频域中看到单个尖峰/峰值的原因。无论正弦波振荡的频率是什么，我们都会在频域中看到尖峰。像这样的尖峰的数学名称称为“脉冲”。

现在如果我们在时域有一个脉冲怎么办？想象一下有人拍手或用锤子敲钉子的录音。这个时频对不太直观。

正如我们所看到的，时域中的尖峰/脉冲在频域中是平坦的，理论上它包含每个频率。理论上不存在完美的脉冲，因为它在时域中必须无限短。与正弦波一样，脉冲发生在时域中的哪个位置并不重要。这里重要的一点是，时域的快速变化会导致出现许多频率。

接下来我们看一下方波的时域和频域图：

这也不太直观，但我们可以看到频域有一个很强的尖峰，恰好位于方波的频率处，但随着频率越高，尖峰就越多。这是由于时域的快速变化，就像前面的例子一样。但它的频率并不平坦。它每隔一段时间就会出现尖峰，并且水平会缓慢衰减（尽管它将永远持续下去）。时域中的方波在频域中具有 sin(x)/x 模式（也称为 sinc 函数）。

现在如果我们在时域中有一个恒定信号怎么办？恒定信号没有“频率”。让我们来看看：

由于没有频率，因此在频域中我们在 0 Hz 处有一个尖峰。如果你仔细想想，这是有道理的。频域不会为“空”，因为这种情况仅在不存在信号时发生（即，时域为 0）。我们将频域中的 0 Hz 称为“DC”，因为它是由及时的 DC 信号（不改变的恒定信号）引起的。请注意，如果我们在时域中增加 DC 信号的幅度，则频域中 0 Hz 处的尖峰也会增加。

稍后我们将了解频域图中 y 轴的确切含义，但现在您可以将其视为一种幅度，告诉您时域信号中存在该频率的多少。

傅里叶变换

从数学上讲，我们用来从时域到频域并返回的“变换”称为傅里叶变换。它的定义如下：

对于信号 x(t)，我们可以使用此公式获得频域版本 X(f)。我们将用 x(t) 或 y(t) 表示函数的时域版本，并用 X(f) 和 Y(f) 表示相应的频域版本。注意“t”代表时间，“f”代表频率。“j”只是虚数单位。您可能在高中数学课上将其视为“i”。我们在工程和计算机科学中使用“j”，因为“i”通常指当前，并且在编程中它通常用作迭代器。

除了比例因子和负号之外，从频率返回到时域几乎是相同的：

请注意，许多教科书和其他资源都使用. 是以弧度每秒为单位的角频率，而是以 Hz 为单位。你所需要知道的是

尽管它在许多方程中添加了一项，但更容易坚持以 Hz 为单位的频率。最终您将在 SDR 应用程序中使用 Hz。

上面的傅立叶变换方程是连续形式，您只会在数学问题中看到它。离散形式更接近于代码中实现的形式：

请注意，主要区别在于我们用求和代替了积分。索引从 0 到 N-1。

如果这些方程对您来说没有多大意义也没关系。实际上，我们不需要直接使用它们来通过 DSP 和 SDR 来做一些很酷的事情！

时频特性

之前我们研究了信号如何在时域和频域中出现的示例。现在，我们将介绍五个重要的“傅立叶性质”。这些属性告诉我们，如果我们对时域信号进行____，那么我们的频域信号就会发生____。它将使我们对在实践中对时域信号执行的数字信号处理 (DSP) 类型有一个重要的了解。

线性特性：

这个属性可能是最容易理解的。如果我们及时将两个信号相加，那么频域版本也将是两个频域信号相加在一起。它还告诉我们，如果我们将其中任何一个乘以缩放因子，则频域也将缩放相同的量。当我们将多个信号加在一起时，此属性的效用将变得更加明显。

频移特性：

x(t) 左边的项就是我们所说的“复正弦曲线”或“复指数”。现在，我们需要知道的是它本质上只是频率为的正弦波。这个属性告诉我们，如果我们获取一个信号并将其乘以正弦波，那么在频域中我们会得到除偏移特定频率之外的。这种频率的转变可能更容易想象：

频移是 DSP 不可或缺的一部分，因为出于多种原因，我们希望上下移动信号的频率。这个属性告诉我们如何做到这一点（乘以正弦波）。这是可视化此属性的另一种方法：

时间缩放属性：

在等式的左侧，我们可以看到我们在时域中缩放信号 x(t)。下面是一个信号随时间缩放的示例，以及每个信号的频域版本会发生什么情况。

时间缩放本质上是缩小或扩大 x 轴上的信号。这个属性告诉我们的是，时域中的缩放会导致频域中的逆缩放。例如，当我们传输比特更快时，我们必须使用更多带宽。该属性有助于解释为什么更高的数据速率信号占用更多的带宽/频谱。如果时间频率缩放是成正比而不是反比，那么蜂窝运营商将能够每秒传输他们想要的所有比特，而无需支付数十亿美元的频谱费用！不幸的是事实并非如此。

那些已经熟悉此属性的人可能会注意到缺少比例因子；为了简单起见，将其省略。出于实际目的，这没有什么区别。

时间属性上的卷积：

之所以称为卷积性质，是因为我们在时域中对 x(t) 和 y(t) 进行卷积。您可能还不了解卷积运算，因此现在将其想象为互相关，尽管我们将在本节中更深入地研究卷积。当我们对时域信号进行卷积时，相当于将这两个信号的频域版本相乘。这与将两个信号加在一起有很大不同。正如我们所见，当您添加两个信号时，实际上什么也没有发生，您只需将频域版本加在一起即可。但是，当您对两个信号进行卷积时，就像从它们创建一个新的第三个信号一样。卷积是 DSP 中最重要的技术，但我们必须首先了解滤波器的工作原理才能完全掌握它。

在我们继续之前，为了简要解释为什么这个属性如此重要，请考虑这种情况：您有一个想要接收的信号，而它旁边有一个干扰信号。

屏蔽的概念在编程中被大量使用，所以我们在这里使用它。如果我们可以创建下面的掩码，并将其乘以上面的信号以掩蔽我们不想要的信号，该怎么办？

我们通常在时域中执行 DSP 操作，因此让我们利用卷积特性来看看如何在时域中进行这种掩蔽。假设 x(t) 是我们接收到的信号。令 Y(f) 为我们要在频域中应用的掩模。这意味着 y(t) 是我们掩模的时域表示，如果我们将它与 x(t) 进行卷积，我们就可以“过滤掉”我们不需要的信号。

当我们讨论过滤时，卷积属性会更有意义。

频率性质的卷积：

最后，我想指出，卷积属性的作用相反，尽管我们不会像时域卷积那样使用它：

还有其他属性，但我认为以上五个属性是最需要理解的。尽管我们没有逐步证明每个属性，但重点是我们使用数学属性来深入了解在进行分析和处理时真实信号会发生什么。不要陷入方程式中。确保您了解每个属性的描述。

快速傅立叶变换 (FFT)

现在回到傅里叶变换。我向您展示了离散傅立叶变换的方程，但您在 99.9% 的编码时间中将使用 FFT 函数 fft()。快速傅里叶变换 (FFT) 只是一种计算离散傅里叶变换的算法。它是几十年前开发的，尽管实现上存在差异，但它仍然是计算离散傅里叶变换的统治者。幸运的是，考虑到他们在名字中使用了“Fast”。

FFT 是一种具有一个输入和一个输出的函数。它将信号从时间转换为频率：

在本教科书中，我们将仅处理一维 FFT（二维用于图像处理和其他应用）。出于我们的目的，可将 FFT 函数视为具有一个输入：样本向量和一个输出：该样本向量的频域版本。输出的大小始终与输入的大小相同。如果我将 1,024 个样本输入 FFT，我将得到 1,024 个样本。令人困惑的部分是，输出将始终位于频域中，因此，如果我们要绘制它，则 x 轴的“跨度”不会根据时域输入中的样本数量而变化。让我们通过查看输入和输出数组及其索引单位来可视化这一点：

由于输出位于频域，因此 x 轴的跨度基于采样率，我们将在下一章介绍。当我们对输入向量使用更多样本时，我们在频域中获得更好的分辨率（除了一次处理更多样本之外）。我们实际上并没有通过更大的输入来“看到”更多的频率。唯一的方法是增加采样率（减少采样周期）。

我们如何实际绘制这个输出？举个例子，假设我们的采样率为每秒 100 万个样本 (1 MHz)。正如我们将在下一章中了解到的，这意味着无论我们向 FFT 提供多少样本，我们只能看到高达 0.5 MHz 的信号。FFT 输出的表示方式如下：

情况总是如此；FFT 的输出将始终显示采样率。即，输出将始终具有负部分和正部分。如果输入很复杂，则负部分和正部分将不同，但如果输入是实数，则它们将相同。

关于频率间隔，每个 bin 对应于Hz，即，向每个 FFT 提供更多样本将导致输出中具有更精细的分辨率。如果您是新手，可以忽略一个非常小的细节：从数学上讲，最后一个索引并不完全对应于，而是对于一个大的值来说大约是。

负频率

负频率到底是什么？现在，只要知道它们与使用复数（虚数）有关即可——在发送/接收射频信号时，实际上并不存在“负频率”之类的东西，它只是我们使用的一种表示形式。这是一种直观的思考方式。假设我们告诉 SDR 调谐到 100 MHz（FM 无线电频段）并以 10 MHz 的速率进行采样。换句话说，我们将查看从 95 MHz 到 105 MHz 的频谱。也许存在三个信号：

现在，当 SDR 给我们样本时，它将如下所示：

请记住，我们将 SDR 调整为 100 MHz。因此，当我们以数字方式表示时，大约 97.5 MHz 的信号显示为 -2.5 MHz，从技术上讲，这是一个负频率。实际上它只是低于中心频率的频率。当我们更多地了解采样并获得使用 SDR 的经验时，这将更有意义。

时间顺序并不重要

在我们开始讨论 FFT 之前，还有最后一个性质。FFT 函数有点“混合”输入信号以形成具有不同比例和单位的输出。我们毕竟不再处于时域中。内化域之间差异的一个好方法是认识到改变时域中发生的事情的顺序不会改变信号中的频率分量。即，以下两个信号的 FFT 都将具有相同的两个尖峰，因为该信号只是不同频率的两个正弦波。改变正弦波发生的顺序并不会改变它们是不同频率的两个正弦波的事实。

从技术上讲，FFT 值的相位会因正弦曲线的时移而发生变化。然而，在本教科书的前几章中，我们将主要关注 FFT 的幅度。

Python 中的 FFT

现在我们已经了解了 FFT 是什么以及如何表示输出，让我们实际查看一些 Python 代码并使用 Numpy 的 FFT 函数 np.fft.fft()。建议您在计算机上使用完整的 Python 控制台/IDE，但在紧要关头，您可以使用左侧导航栏底部链接的基于 Web 的在线 Python 控制台。

首先我们需要在时域中创建一个信号。请随意使用您自己的 Python 控制台。为了简单起见，我们将制作一个 0.15 Hz 的简单正弦波。我们还将使用 1 Hz 的采样率，这意味着我们在 0、1、2、3 秒等时间进行采样。

import numpy as np
t = np.arange(100)
s = np.sin(0.15*2*np.pi*t)

如果我们绘制s它看起来像：

接下来我们使用 Numpy 的 FFT 函数：

S = np.fft.fft(s)

如果我们看一下，S我们会发现它是一个复数数组：

S =  array([-0.01865008 +0.00000000e+00j, -0.01171553 -2.79073782e-01j,0.02526446 -8.82681208e-01j,  3.50536075 -4.71354150e+01j, -0.15045671 +1.31884375e+00j, -0.10769903 +7.10452463e-01j, -0.09435855 +5.01303240e-01j, -0.08808671 +3.92187956e-01j, -0.08454414 +3.23828386e-01j, -0.08231753 +2.76337148e-01j, -0.08081535 +2.41078885e-01j, -0.07974909 +2.13663710e-01j,...

提示：无论您在做什么，如果遇到复数，请尝试计算幅度和相位，看看它们是否更有意义。让我们准确地做到这一点，并绘制幅度和相位。在大多数语言中，abs() 是表示复数大小的函数。相位的函数各不相同，但在 Python 中是np.angle().

import matplotlib.pyplot as plt
S_mag = np.abs(S)
S_phase = np.angle(S)
plt.plot(t,S_mag,'.-')
plt.plot(t,S_phase,'.-')

现在我们没有为绘图提供任何 x 轴，它只是数组的索引（从 0 开始计数）。由于数学原因，FFT 的输出具有以下格式：

但我们希望 0 Hz (DC) 位于中心，负频率位于左侧（这正是我们希望将事物可视化的方式）。因此，每当我们进行 FFT 时，我们都需要执行“FFT 移位”，这只是一个简单的数组重新排列操作，有点像循环移位，但更多的是“把这个放在那里”。下图完整定义了 FFT 移位运算的作用：

为了方便起见，Numpy 有一个 FFT 移位函数np.fft.fftshift()。将 np.fft.fft() 行替换为：

S = np.fft.fftshift(np.fft.fft(s))

我们还需要计算出 x 轴值/标签。回想一下，为了简单起见，我们使用了 1 Hz 的采样率。这意味着频域图的左边缘将为 -0.5 Hz，右边缘将为 0.5 Hz。如果这没有意义，那么在您读完IQ 采样一章后就会明白了。让我们坚持采样率为 1 Hz 的假设，并使用适当的 x 轴标签绘制 FFT 输出的幅度和相位。以下是此 Python 示例的最终版本和输出：

import numpy as np
import matplotlib.pyplot as plt

Fs = 1 # Hz
N = 100 # number of points to simulate, and our FFT size

t = np.arange(N) # because our sample rate is 1 Hz
s = np.sin(0.15*2*np.pi*t)
S = np.fft.fftshift(np.fft.fft(s))
S_mag = np.abs(S)
S_phase = np.angle(S)
f = np.arange(Fs/-2, Fs/2, Fs/N)
plt.figure(0)
plt.plot(f, S_mag,'.-')
plt.figure(1)
plt.plot(f, S_phase,'.-')
plt.show()

请注意，我们看到峰值为 0.15 Hz，这是我们创建正弦波时使用的频率。这意味着我们的 FFT 成功了！如果我们不知道用于生成该正弦波的代码，但只给出了样本列表，则可以使用 FFT 来确定频率。我们之所以在 -0.15 Hz 处看到尖峰，是因为它是一个真实的信号，而不是复杂的信号，稍后我们将更深入地讨论这一点。

窗口化

当我们使用 FFT 测量信号的频率分量时，FFT 假设给它的是一段周期信号。它的行为就好像我们提供的信号片段继续无限重复。就好像切片的最后一个样本连接回第一个样本。它源于傅里叶变换背后的理论。这意味着我们希望避免第一个和最后一个样本之间的突然转变，因为时域中的突然转变看起来像许多频率，而实际上我们的最后一个样本实际上并没有连接回我们的第一个样本。简而言之：如果我们要使用进行 100 个样本的 FFT，np.fft.fft(x)我们希望x[0]和x[99]的值相等或接近。

我们弥补这种循环特性的方法是通过“窗口”。在 FFT 之前，我们将信号切片乘以窗函数，该函数就是两端逐渐变细到零的函数。这确保了信号片段将以零开始和结束并连接。常见的窗函数包括 Hamming、Hanning、Blackman 和 Kaiser。当您不应用任何窗口时，称为使用“矩形”窗口，因为它就像乘以一个数组。以下是几个窗口函数的样子：

对于初学者来说，一个简单的方法是坚持使用汉明窗，它可以在 Python 中创建，其中np.hamming(N)N 是数组中元素的数量，即 FFT 大小。在上面的练习中，我们将在 FFT 之前应用窗口。在第二行代码之后我们将插入：

s = s * np.hamming(100)

如果您担心选择错误的窗户，请不要担心。与根本不使用窗口相比，Hamming、Hanning、Blackman 和 Kaiser 之间的差异非常小，因为它们在两侧都逐渐变细到零并解决了根本问题。

FFT 大小调整

最后要注意的是 FFT 大小调整。由于 FFT 的实现方式，最佳 FFT 大小始终为 2 的量级。您可以使用不是 2 数量级的大小，但速度会更慢。常见尺寸在 128 到 4,096 之间，当然您也可以选择更大的尺寸。在实践中，我们可能需要处理数百万或数十亿个样本长的信号，因此我们需要分解信号并进行多次 FFT。这意味着我们将得到很多输出。我们可以将它们平均或随时间绘制它们（特别是当我们的信号随时间变化时）。您不必将信号的每个样本都通过 FFT 来获得该信号的良好频域表示。例如，您只能对信号中的每 100k 个样本进行 1,024 次 FFT，但只要信号始终开启，它可能看起来仍然不错。

频谱图/瀑布图

频谱图是显示频率随时间变化的图。它只是一堆堆叠在一起的 FFT（如果您想要水平轴上的频率，则垂直堆叠）。我们还可以实时展示它，通常称为瀑布。频谱分析仪是显示该频谱图/瀑布图的设备。下图显示了如何将 IQ 样本数组分割以形成频谱图：

由于频谱图涉及绘制 2D 数据，因此它实际上是 3D 绘图，因此我们必须使用颜色图来表示 FFT 幅值，即我们想要绘制的“值”。以下是频谱图的示例，水平/x 轴为频率，垂直/y 轴为时间。蓝色代表最低能量，红色代表最高能量。我们可以看到中心有一个强烈的 DC (0 Hz) 尖峰，周围有变化的信号。蓝色代表我们的本底噪声。

请记住，它只是相互堆叠的 FFT 行，每行是 1 个 FFT（从技术上讲，1 个 FFT 的量级）。请务必将输入信号按 FFT 大小的切片进行时间切片（例如，每切片 1024 个样本）。在进入代码生成频谱图之前，这是我们将使用的示例信号，它只是白噪声中的一个音调：

import numpy as np
import matplotlib.pyplot as plt

sample_rate = 1e6

# Generate tone plus noise
t = np.arange(1024*1000)/sample_rate # time vector
f = 50e3 # freq of tone
x = np.sin(2*np.pi*f*t) + 0.2*np.random.randn(len(t))

以下是时域中的情况（前 200 个样本）：

在 Python 中，我们可以生成频谱图，如下所示：

# simulate the signal above, or use your own signal

fft_size = 1024
num_rows = len(x) // fft_size # // is an integer division which rounds down
spectrogram = np.zeros((num_rows, fft_size))
for i in range(num_rows):
    spectrogram[i,:] = 10*np.log10(np.abs(np.fft.fftshift(np.fft.fft(x[i*fft_size:(i+1)*fft_size])))**2)

plt.imshow(spectrogram, aspect='auto', extent = [sample_rate/-2/1e6, sample_rate/2/1e6, len(x)/sample_rate, 0])
plt.xlabel("Frequency [MHz]")
plt.ylabel("Time [s]")
plt.show()

这应该产生以下结果，这不是最有趣的频谱图，因为没有随时间变化的行为。有两种音调，因为我们模拟了真实信号，而真实信号始终具有与正侧相匹配的负 PSD。有关频谱图的更多有趣示例，请查看https://www.IQEngine.org！

FFT 实现

尽管 NumPy 已经为我们实现了 FFT，但了解其底层工作原理还是很高兴的。最流行的 FFT 算法是 Cooley-Tukey FFT 算法，最初由 Carl Friedrich Gauss 于 1805 年左右发明，后来由 James Cooley 和 John Tukey 于 1965 年重新发现并推广。

该算法的基本版本适用于二次方大小的 FFT，适用于复杂输入，但也适用于实际输入。该算法的构建块被称为蝴蝶，它本质上是一个 N = 2 大小的 FFT，由两个乘法和两个求和组成：

或者

其中被称为旋转因子（是子 FFT 的大小，是索引）。请注意，输入和输出应该是复数，例如，可能是 0.6123 - 0.5213j，并且和/乘法是复数。

该算法是递归的，并将自身分成两半，直到剩下一系列蝴蝶，下面使用大小为 8 的 FFT 对此进行了描述：

此模式中的每一列都是一组可以并行完成的操作，并且执行步骤，这就是为什么 FFT 的计算复杂度是 O(· )，而 DFT 是 O(· )。

对于那些更喜欢用代码而不是方程式思考的人，下面显示了 FFT 的简单 Python 实现，以及一个由音调和噪声组成的示例信号，用于尝试 FFT。

import numpy as np
import matplotlib.pyplot as plt

def fft(x):
    N = len(x)
    if N == 1:
        return x
    twiddle_factors = np.exp(-2j * np.pi * np.arange(N//2) / N)
    x_even = fft(x[::2]) # yay recursion!
    x_odd = fft(x[1::2])
    return np.concatenate([x_even + twiddle_factors * x_odd,
                           x_even - twiddle_factors * x_odd])

# Simulate a tone + noise
sample_rate = 1e6
f_offset = 0.2e6 # 200 kHz offset from carrier
N = 1024
t = np.arange(N)/sample_rate
s = np.exp(2j*np.pi*f_offset*t)
n = (np.random.randn(N) + 1j*np.random.randn(N))/np.sqrt(2) # unity complex noise
r = s + n # 0 dB SNR

# Perform fft, fftshift, convert to dB
X = fft(r)
X_shifted = np.roll(X, N//2) # equivalent to np.fft.fftshift
X_mag = 10*np.log10(np.abs(X_shifted)**2)

# Plot results
f = np.linspace(sample_rate/-2, sample_rate/2, N)/1e6 # plt in MHz
plt.plot(f, X_mag)
plt.plot(f[np.argmax(X_mag)], np.max(X_mag), 'rx') # show max
plt.grid()
plt.xlabel('Frequency [MHz]')
plt.ylabel('Magnitude [dB]')
plt.show()