DSP - 如何处理 IQ 数据
前言
在上一篇文章中，我们介绍了解析 IQ 数据的方法与思路。本文将探讨 IQ 数据的 FFT (快速傅里叶变换) 以及频域上的截取操作，阅读本文可以学习到以下内容：
- 如何将 IQ 数据通过 FFT 转换为频域数据。
- 如何实现 IQ 文件的频域截取。
- 如何将截取后的频域数据转回 IQ 数据 (时域)。
- 如何将复数转回实数对。
将 IQ 数据通过 FFT 转换为频域数据
我们将实现解析并处理一个以 int16 实数对存储的 IQ 文件，并对该文件的 IQ 数据进行 FFT 转换为频域数据，最终绘制一个频谱图展示该信号。
import os
import numpy as np
import matplotlib.pyplot as plt
def read_iq_file_in_chunks(filename, dtype=np.int16, chunk_size=2048 * 200):
sample_size = np.dtype(dtype).itemsize
read_len = chunk_size * sample_size * 2
with open(filename, 'rb') as f:
while True:
data = f.read(read_len) # 每个样本占 4 字节 (2 字节 I + 2 字节 Q)
if not data:
break
iq_data = np.frombuffer(data, dtype=np.int16)
yield iq_data
def compute_spectrum(iq_data, fs, fft_size=2048):
samples = iq_data[::2] + 1j * iq_data[1::2]
segments = range(len(samples) // fft_size)
for i in segments:
fft_window = np.hanning(fft_size)
windowed_data = samples[fft_size * i:fft_size * (i + 1)] * fft_window
# 做完 fft 后将 oHz 频率挪至中心
fft_result = np.fft.fftshift(np.fft.fft(windowed_data, n=fft_size))
fft_freq = np.fft.fftshift(np.fft.fftfreq(fft_size, 1 / fs))
# 计算功率
magnitude_spectrum = np.abs(fft_result)
PSD = magnitude_spectrum ** 2 / (fft_size * fs)
PSD_log = 10.0 * np.log10(PSD)
# 绘制频谱图
plt.plot(fft_freq, PSD_log)
plt.xlabel("Frequency [Hz]")
plt.ylabel("Magnitude [dB]")
plt.grid(True)
plt.show()
def main():
fft_size = 2048
dtype = np.int16
file = r'F:\SignalFiles\1000000000_800000_1280000.iq'
file_name = os.path.basename(file)
names = os.path.splitext(file_name)[0].split('_')
center_freq = float(names[0])
bandwidth = float(names[1])
sample_rate = float(names[2])
for iq_data in read_iq_file_in_chunks(file, dtype):
compute_spectrum(iq_data, sample_rate, fft_size)
if __name__ == "__main__":
main()
接下来我们解释一下这段代码都做了哪些事情:
- 首先定义了分段读取 IQ 文件的函数
read_iq_file_in_chunks()和计算频谱的函数compute_spectrum - 在
main()方法中，定义了我们做 FFT 的窗大小为 2048，IQ 文件的存储格式为int16，然后解析这个 IQ 文件名得到中心频率、带宽和采样率参数。 - 接下来我们分段读取 IQ 数据，由于我设置的 FFT 窗大小为 2048，所以我做一次 FFT 需要使用 2048 个采样点对应频谱图的一帧数据，我将一次性取出 200 帧采样点进行处理。
- 然后我们将读取出来的采样点数据、采样率以及FFT 窗大小参数传入
compute_spectrum()函数中。 - 在
compute_spectrum()函数中，我们先将实数转为复数，然后一帧一帧的进行 FFT 转换。 - FFT 窗函数使用的是
Hanning窗，使用np.fft.fft()将 IQ 数据转换为频谱数据，然后再使用np.fft.fftshift()将负频率部分搬移至最左侧。同时使用np.fft.fftfreq()获取频率轴。 - 为了将信号展示在频谱图上，我们需要计算功率并将单位转换为 dB，最后将数据绘制出来得到如下图所示的图像：
实现 IQ 文件的频域截取
为什么要做频域截取
在信号分析业务中，我们经常要对检测到的信号做截取，或者对某一频段做截取。用上面这张图举例，由于该信号文件的采样率为 1.28 MHz/s, 而采样带宽为 800kHz ，这就导致在频域上两边会有一个下沿，而我们需要的数据实际上是采样带宽范围内的 IQ 数据！如下图所示：
频域截取
接下来我们进行频域的截取流程：
import os
import numpy as np
import matplotlib.pyplot as plt
def read_iq_file_in_chunks(filename, dtype=np.int16, chunk_size=2048 * 200):
sample_size = np.dtype(dtype).itemsize
read_len = chunk_size * sample_size * 2
with open(filename, 'rb') as f:
while True:
data = f.read(read_len) # 每个样本占 4 字节 (2 字节 I + 2 字节 Q)
if not data:
break
iq_data = np.frombuffer(data, dtype=np.int16)
yield iq_data
def compute_spectrum(iq_data, fs, bandwidth, fft_size=2048):
samples = iq_data[::2] + 1j * iq_data[1::2]
segments = range(len(samples) // fft_size)
for i in segments:
fft_window = np.hanning(fft_size)
windowed_data = samples[fft_size * i:fft_size * (i + 1)] * fft_window
# 做完 fft 后将 oHz 频率挪至中心
fft_result = np.fft.fftshift(np.fft.fft(windowed_data, n=fft_size))
fft_freq = np.fft.fftshift(np.fft.fftfreq(fft_size, 1 / fs))
# ----- 新增 ----
# 截取真实带宽频域数据
center_idx = fft_size // 2
half_bw_idx = int(bandwidth / 2 * fft_size / fs)
start_idx = center_idx - half_bw_idx
end_idx = center_idx + half_bw_idx
# 频域截取
signal_fft = np.zeros_like(fft_result)
signal_fft[start_idx:end_idx] = fft_result[start_idx:end_idx]
# --------------
# 将复数转为实数
magnitude_spectrum = np.abs(signal_fft)
# 计算功率
PSD = magnitude_spectrum ** 2 / (fft_size * fs)
PSD_log = 10.0 * np.log10(PSD)
# 画图
plt.plot(fft_freq, PSD_log)
plt.xlabel("Frequency [Hz]")
plt.ylabel("Magnitude [dB]")
plt.grid(True)
plt.show()
print("success")
def main():
fft_size = 2048
dtype = np.int16
file = r'F:\SignalFiles\1000000000_800000_1280000.iq'
file_name = os.path.basename(file)
names = os.path.splitext(file_name)[0].split('_')
center_freq = float(names[0])
bandwidth = float(names[1])
sample_rate = float(names[2])
for iq_data in read_iq_file_in_chunks(file, dtype):
compute_spectrum(iq_data, sample_rate, bandwidth, fft_size)
if __name__ == "__main__":
main()
接下来我们解释一下这段代码的运行流程，大部分和 FFT 转换那部分的代码一致，只是新增了截取的操作。
- 在进行 FFT 以后，由于我们将负频率部分搬移到了最左侧，所以中心频率所在的索引点就是该数组长度的一半也就是
fft_size / 2。 - 接下来我们算出带宽一半的索引值，再算出我们需要数据的索引起始点和结束点。
- 最后我们先创建一个长度一致的都为 0 的数组，再将索引起始点至结束点的值拷贝过来。
- 至此，我们就成功的进行了频域的截取，得到了我们想要的数据，如下图所示：
实现频域数据转回 IQ 数据 (时域)
当我们从频域截取完成后，还需要再将频域数据逆 FFT 转回时域，因为我们可能需要用到截取后的 IQ 数据，比如进行信号识别等业务。 实现方法非常简单：
# 转换回时域
# signal_fft 是频域截取后的数据
extracted_signal_time = np.fft.ifft(np.fft.ifftshift(signal_fft))
real_part = np.real(extracted_signal_time)
imag_part = np.imag(extracted_signal_time)
# 交错合并为实数数组
iq_real = np.empty(extracted_signal_time.size * 2, dtype=np.float32)
iq_real[0::2] = real_part
iq_real[1::2] = imag_part
iq_data = np.real(iq_real) # 实数 IQ 数据
- 首先我们需要将截取后的频域数据使用
np.fft.ifftshift()搬回之前的位置，还记得吗，我们之前把负频率的部分搬到了最左侧，现在需要搬回去。 - 使用
np.fft.ifft()逆向 FFT 将频域数据转回 IQ 数据。 - 提取复数的实部和虚部然后交错合并为实数数组。
总结
在本文中，我们详细探讨了如何将IQ数据通过快速傅里叶变换（FFT）转换为频域数据。我们首先介绍了IQ文件的频域提取的必要性，分析了提取频域数据所带来的优势和应用场景。
接着，我们深入讨论了频域数据提取的具体步骤，提供了实用的示例和技巧。最后，我们探讨了如何将频域数据转换回IQ数据（时域），并将复数转换为实数对。
尽管我在本文中尽可能地详细介绍了每一个步骤和细节，但是难免会存在一些错误和不足之处。如果您在使用本文中介绍的方法时发现了任何错误或者有更好的方法，非常欢迎您指正并提出建议，以便我能够不断改进和提升文章的质量。
我是Ricky，一个兴趣使然的开发者。非常感谢您阅读本文，希望本文对您有所帮助！