Hướng dẫn python get frequency of signal - python lấy tần số tín hiệu

Không rõ câu hỏi của bạn chính xác những giá trị trong tệp đại diện. Nhưng giả sử rằng chúng chỉ ra các mẫu điện áp liên tiếp, bạn có thể tải tệp vào một mảng numpy bằng cách sử dụng

import numpy as np
data = np.array([float(f) for f in file(filename).read().split()])

và sau đó tính toán biến đổi Fourier là

import numpy.fft as fft
spectrum = fft.fft(data)

Sau đó, bạn có thể vẽ đồ thị cường độ của FFT là

freq = fft.fftfreq(len(spectrum))
plot(freq, abs(spectrum))

Và những gì bạn thấy nên phù hợp với những gì hiển thị trên máy hiện sóng.

Nếu bạn muốn xác định các tần số chiếm ưu thế trong phổ, bạn sẽ phải cắt mảng ở một số ngưỡng, ví dụ: một cái gì đó như thế này:

threshold = 0.5 * max(abs(spectrum))
mask = abs(spectrum) > threshold
peaks = freq[mask]

Nội dung của

# Make plots appear inline, set custom plotting style
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('style/elegant.mplstyle')

4 (và do đó cũng

# Make plots appear inline, set custom plotting style
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('style/elegant.mplstyle')

5) là tần số trong các đơn vị của tốc độ lấy mẫu. Ví dụ: nếu máy hiện sóng của bạn lấy mẫu dạng sóng mỗi micro giây, các giá trị trong

# Make plots appear inline, set custom plotting style
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('style/elegant.mplstyle')

4 nằm trong megahertz. Vì vậy, nếu bạn cung cấp tín hiệu 1 kHz lý tưởng, bạn có thể làm với ví dụ:

t = arange(4e6) / 1e6 # sampling times in seconds
data = sin(2 * pi * 1000 * t)

Bạn sẽ nhận được một đỉnh cao ở mức 0,001 MHz, và theo đó bạn sẽ thấy

# Make plots appear inline, set custom plotting style
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('style/elegant.mplstyle')

7.

Chương 4. Tần suất và biến đổi Fourier nhanhFrequency and the Fast Fourier Transform

Nếu bạn muốn tìm những bí mật của vũ trụ, hãy suy nghĩ về năng lượng, tần suất và rung động.

Nikola Tesla

Chương này được viết với sự hợp tác của cha SW, PW Van der Walt.

Chương này sẽ khởi hành một chút từ định dạng của phần còn lại của cuốn sách. Cụ thể, bạn có thể tìm thấy mã trong chương khá khiêm tốn. Thay vào đó, chúng tôi muốn minh họa một thuật toán thanh lịch, biến đổi Fourier nhanh (FFT), rất hữu ích, được thực hiện trong SCIPY, và tất nhiên, trên các mảng Numpy.

Giới thiệu tần số

Chúng tôi sẽ bắt đầu bằng cách thiết lập một số kiểu vẽ và nhập các nghi phạm thông thường:

# Make plots appear inline, set custom plotting style
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('style/elegant.mplstyle')

import numpy as np

Biến đổi Fourier rời rạc (DFT) là một kỹ thuật toán học được sử dụng để chuyển đổi dữ liệu thời gian hoặc không gian thành dữ liệu miền tần số. Tần số là một khái niệm quen thuộc, do sự xuất hiện thông tục của nó trong ngôn ngữ tiếng Anh: Ghi chú thấp nhất mà tai nghe của bạn có thể ầm ầm là khoảng 20 Hz, trong khi C giữa trên cây đàn piano nằm khoảng 261,6 Hz; Hertz, hoặc dao động mỗi giây, trong trường hợp này thực sự đề cập đến số lần mỗi giây mà màng bên trong tai nghe di chuyển đến và-fro. Điều đó, đến lượt nó, tạo ra các xung không khí nén, khi đến màng nhĩ của bạn, gây ra rung động ở cùng tần số. Vì vậy, nếu bạn thực hiện một chức năng định kỳ đơn giản, sin (10 × 2πt), bạn có thể xem nó dưới dạng sóng:

f = 10  # Frequency, in cycles per second, or Hertz
f_s = 100  # Sampling rate, or number of measurements per second

t = np.linspace(0, 2, 2 * f_s, endpoint=False)
x = np.sin(f * 2 * np.pi * t)

fig, ax = plt.subplots()
ax.plot(t, x)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Signal amplitude');

Hoặc bạn có thể nghĩ tương đương với nó như một tín hiệu lặp lại tần số 10 Hz (nó lặp lại một lần 1/10 giây một khoảng thời gian chúng ta gọi là kỳ của nó). Mặc dù chúng tôi tự nhiên liên kết tần số theo thời gian, nhưng nó cũng có thể được áp dụng cho không gian. Ví dụ, một bức ảnh của một mẫu dệt thể hiện tần số không gian cao, trong khi bầu trời hoặc các vật thể mịn khác có tần số không gian thấp.

Bây giờ chúng ta hãy kiểm tra xoang của chúng ta thông qua ứng dụng DFT:

from scipy import fftpack

X = fftpack.fft(x)
freqs = fftpack.fftfreq(len(x)) * f_s

fig, ax = plt.subplots()

ax.stem(freqs, np.abs(X))
ax.set_xlabel('Frequency in Hertz [Hz]')
ax.set_ylabel('Frequency Domain (Spectrum) Magnitude')
ax.set_xlim(-f_s / 2, f_s / 2)
ax.set_ylim(-5, 110)

(-5, 110)

Chúng ta thấy rằng đầu ra của FFT là một mảng 1D có hình dạng giống như đầu vào, chứa các giá trị phức tạp. Tất cả các giá trị đều bằng không, ngoại trừ hai mục. Theo truyền thống, chúng tôi hình dung độ lớn của kết quả như một biểu đồ gốc, trong đó chiều cao của mỗi thân tương ứng với giá trị cơ bản.

(Chúng tôi giải thích lý do tại sao bạn thấy các tần số tích cực và tiêu cực sau này trong các biến đổi Fourier rời rạc.

Biến đổi Fourier đưa chúng ta từ thời điểm đến miền tần số và điều này hóa ra là có một số lượng lớn các ứng dụng. Biến đổi Fourier nhanh (FFT) là một thuật toán để tính toán DFT; Nó đạt được tốc độ cao của nó bằng cách lưu trữ và tái sử dụng kết quả tính toán khi nó tiến triển.

Trong chương này, chúng tôi kiểm tra một vài ứng dụng của DFT để chứng minh rằng FFT có thể được áp dụng cho dữ liệu đa chiều (không chỉ các phép đo 1D) để đạt được nhiều mục tiêu.

Minh họa: Một phổ chim Birdsong

Hãy bắt đầu với một trong những ứng dụng phổ biến nhất, chuyển đổi tín hiệu âm thanh (bao gồm các biến thể của áp suất không khí theo thời gian) thành một phổ. Bạn có thể đã thấy các quang phổ trên chế độ xem bộ cân bằng trình phát nhạc của bạn, hoặc thậm chí trên một âm thanh nổi trường học cũ (Hình & NBSP; 4-1).

Nghe đoạn trích này của The Nightingale Birdsong (phát hành theo CC vào 4.0):

import numpy.fft as fft
spectrum = fft.fft(data)

0

Nếu bạn đang đọc phiên bản giấy của cuốn sách này, bạn sẽ phải sử dụng trí tưởng tượng của mình! Nó đi một cái gì đó như thế này: Chee-Chee-Woorrrr-hee-hee Cheet-Wheet-Hoorrr-Chirr-wiki-wheo-wheo-wheo-wheo-wheo-wheo-wheo.

Vì chúng tôi nhận ra rằng không phải ai cũng thông thạo tiếng chim, có lẽ nó tốt nhất nếu chúng tôi hình dung ra các phép đo mà Better được gọi là tín hiệu của Hồi giáo.

Chúng tôi tải tệp âm thanh, cung cấp cho chúng tôi tốc độ lấy mẫu (số lượng phép đo mỗi giây) cũng như dữ liệu âm thanh dưới dạng mảng

# Make plots appear inline, set custom plotting style
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('style/elegant.mplstyle')

8 hai cột vì đây là bản ghi âm thanh nổi.

import numpy.fft as fft
spectrum = fft.fft(data)

1

Chúng tôi chuyển đổi thành Mono bằng cách lấy trung bình các kênh bên trái và bên phải.

import numpy.fft as fft
spectrum = fft.fft(data)

2

Sau đó, chúng tôi tính toán độ dài của đoạn trích và vẽ âm thanh (Hình & NBSP; 4-2).

import numpy.fft as fft
spectrum = fft.fft(data)

3

import numpy.fft as fft
spectrum = fft.fft(data)

4

Hình 4-2. Biểu đồ dạng sóng âm thanh của Nightingale BirdsongAudio waveform plot of nightingale birdsong

Chà, đó không phải là rất thỏa mãn, phải không? Nếu tôi gửi điện áp này cho một loa, tôi có thể nghe thấy một con chim kêu vang, nhưng tôi có thể tưởng tượng rất rõ nó sẽ phát ra âm thanh trong đầu tôi. Có cách nào tốt hơn để xem những gì đang xảy ra?

Có, và nó được gọi là biến đổi Fourier rời rạc, hoặc DFT, trong đó rời rạc đề cập đến bản ghi bao gồm các phép đo âm thanh cách nhau thời gian, trái ngược với bản ghi liên tục như, ví dụ, trên băng từ tính (có nhớ băng cassette không?). DFT thường được tính toán bằng thuật toán FFT, một tên được sử dụng không chính thức để chỉ chính DFT. DFT cho chúng ta biết tần số hoặc ghi chú của người dùng để mong đợi trong tín hiệu của chúng ta.

Tất nhiên, một con chim hát nhiều nốt nhạc trong suốt bài hát, vì vậy chúng tôi cũng muốn biết khi nào mỗi nốt xảy ra. Biến đổi Fourier lấy tín hiệu trong miền thời gian (nghĩa là, một tập hợp các phép đo theo thời gian) và biến nó thành một phổ, một tập hợp các tần số có các giá trị tương ứng (phức tạp2). Phổ không chứa bất kỳ thông tin nào về thời gian! 3

Vì vậy, để tìm ra cả tần số và thời gian mà chúng được hát, chúng tôi sẽ cần phải có phần thông minh. Chiến lược của chúng tôi như sau: Lấy tín hiệu âm thanh, chia nó thành các lát nhỏ, chồng chéo và áp dụng biến đổi Fourier cho mỗi (một kỹ thuật được gọi là biến đổi Fourier thời gian ngắn).

Chúng tôi sẽ chia tín hiệu thành các lát 1.024 mẫu, đó là khoảng 0,02 giây âm thanh. Lý do chúng tôi đã chọn 1.024 và không phải 1.000 chúng tôi sẽ giải thích trong một giây khi chúng tôi kiểm tra hiệu suất. Các lát cắt sẽ chồng chéo bởi 100 mẫu như được hiển thị ở đây:

Bắt đầu bằng cách cắt tín hiệu thành các lát của 1024 mẫu, mỗi lát chồng lên nhau bằng 100 mẫu trước đó. Đối tượng

# Make plots appear inline, set custom plotting style
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('style/elegant.mplstyle')

9 kết quả chứa một lát mỗi hàng.

import numpy.fft as fft
spectrum = fft.fft(data)

5

import numpy.fft as fft
spectrum = fft.fft(data)

6

Tạo một chức năng cửa sổ (xem phần Winding Winding để thảo luận về các giả định và diễn giải cơ bản của từng người) và nhân nó với tín hiệu:

import numpy.fft as fft
spectrum = fft.fft(data)

7

Nó thuận tiện hơn để có một lát mỗi cột, vì vậy chúng tôi thực hiện chuyển vị:

import numpy.fft as fft
spectrum = fft.fft(data)

8

import numpy.fft as fft
spectrum = fft.fft(data)

9

Đối với mỗi lát, tính toán DFT, trả về cả tần số dương và âm (nhiều hơn về tần số của các tần số và thứ tự của chúng), vì vậy chúng tôi sẽ cắt các tần số M2 dương.

freq = fft.fftfreq(len(spectrum))
plot(freq, abs(spectrum))

0

.

Phổ có thể chứa cả giá trị rất lớn và rất nhỏ. Lấy nhật ký nén phạm vi đáng kể.

Ở đây chúng tôi thực hiện một biểu đồ nhật ký của tỷ lệ tín hiệu chia cho tín hiệu tối đa (được hiển thị trong Hình & NBSP; 4-3). Đơn vị cụ thể được sử dụng cho tỷ lệ là decibel, 20Log10 (tỷ lệ biên độ).

freq = fft.fftfreq(len(spectrum))
plot(freq, abs(spectrum))

1

Hình 4-3. Phổ chimBirdsong spectrogram

Tốt hơn nhiều! Bây giờ chúng ta có thể thấy rằng các tần số thay đổi theo thời gian và quang phổ tương ứng với cách âm thanh âm thanh. Xem nếu bạn có thể phù hợp với mô tả trước đây của chúng tôi: Chee-Chee-Woorrrrrr-hee-hee Cheet-Wheet-Hoorrrrr-wiki-wheo-wheo-wheo-wheo-wheo-wheo-wheo. .

SCIPY đã bao gồm việc thực hiện quy trình này là

import numpy as np

3 (Hình & NBSP; 4-4), có thể được gọi như sau:

freq = fft.fftfreq(len(spectrum))
plot(freq, abs(spectrum))

2

Hình 4-4. SCIPY TÌM KIẾM TIÊU CHUẨN CỦA BIRDSONG CHUYỂN ĐỔISciPy built-in rendition of birdsong spectrogram

Sự khác biệt duy nhất giữa quang phổ thủ công mà chúng tôi đã tạo so với hàm tích hợp SCIPY là SCIPY trả về bình phương cường độ phổ (biến điện áp đo thành năng lượng đo được) và nhân lên nó với một số yếu tố bình thường hóa.4

Lịch sử

Theo dõi nguồn gốc chính xác của biến đổi Fourier là khó khăn. Một số thủ tục liên quan đi xa như thời Babylon, nhưng đó là các chủ đề nóng của việc tính toán quỹ đạo tiểu hành tinh và giải phương trình nhiệt (dòng chảy) dẫn đến một số đột phá vào đầu những năm 1800. Ai chính xác giữa Clairaut, LaGrange, Euler, Gauss và D'Alembert Chúng ta nên cảm ơn không chính xác, nhưng Gauss là người đầu tiên mô tả biến đổi Fourier nhanh (một thuật toán để tính toán DFT, được phổ biến bởi Cooley và Tukey vào năm 1965) . Joseph Fourier, sau khi biến đổi được đặt tên, lần đầu tiên tuyên bố rằng các hàm định kỳ5 tùy ý có thể được biểu thị bằng tổng số các hàm lượng giác.

Thực hiện

Chức năng DFT trong SCIPY sống trong mô -đun

import numpy as np

4. Trong số những thứ khác, nó cung cấp chức năng liên quan đến DFT sau:

import numpy as np

2,

import numpy as np

6,

import numpy as np

7

Tính toán DFT bằng thuật toán FFT trong kích thước 1, 2 hoặc

import numpy as np

8.

import numpy as np

9,

f = 10  # Frequency, in cycles per second, or Hertz
f_s = 100  # Sampling rate, or number of measurements per second

t = np.linspace(0, 2, 2 * f_s, endpoint=False)
x = np.sin(f * 2 * np.pi * t)

fig, ax = plt.subplots()
ax.plot(t, x)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Signal amplitude');

0,

f = 10  # Frequency, in cycles per second, or Hertz
f_s = 100  # Sampling rate, or number of measurements per second

t = np.linspace(0, 2, 2 * f_s, endpoint=False)
x = np.sin(f * 2 * np.pi * t)

fig, ax = plt.subplots()
ax.plot(t, x)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Signal amplitude');

1

Tính nghịch đảo của DFT.

f = 10  # Frequency, in cycles per second, or Hertz
f_s = 100  # Sampling rate, or number of measurements per second

t = np.linspace(0, 2, 2 * f_s, endpoint=False)
x = np.sin(f * 2 * np.pi * t)

fig, ax = plt.subplots()
ax.plot(t, x)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Signal amplitude');

2,

f = 10  # Frequency, in cycles per second, or Hertz
f_s = 100  # Sampling rate, or number of measurements per second

t = np.linspace(0, 2, 2 * f_s, endpoint=False)
x = np.sin(f * 2 * np.pi * t)

fig, ax = plt.subplots()
ax.plot(t, x)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Signal amplitude');

3,

f = 10  # Frequency, in cycles per second, or Hertz
f_s = 100  # Sampling rate, or number of measurements per second

t = np.linspace(0, 2, 2 * f_s, endpoint=False)
x = np.sin(f * 2 * np.pi * t)

fig, ax = plt.subplots()
ax.plot(t, x)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Signal amplitude');

4,

f = 10  # Frequency, in cycles per second, or Hertz
f_s = 100  # Sampling rate, or number of measurements per second

t = np.linspace(0, 2, 2 * f_s, endpoint=False)
x = np.sin(f * 2 * np.pi * t)

fig, ax = plt.subplots()
ax.plot(t, x)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Signal amplitude');

5

Tính toán các biến đổi cosin và sin, và nghịch đảo của chúng.

f = 10  # Frequency, in cycles per second, or Hertz
f_s = 100  # Sampling rate, or number of measurements per second

t = np.linspace(0, 2, 2 * f_s, endpoint=False)
x = np.sin(f * 2 * np.pi * t)

fig, ax = plt.subplots()
ax.plot(t, x)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Signal amplitude');

6,

f = 10  # Frequency, in cycles per second, or Hertz
f_s = 100  # Sampling rate, or number of measurements per second

t = np.linspace(0, 2, 2 * f_s, endpoint=False)
x = np.sin(f * 2 * np.pi * t)

fig, ax = plt.subplots()
ax.plot(t, x)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Signal amplitude');

7

Chuyển thành phần tần số bằng không thành trung tâm của phổ và quay lại, tương ứng (nhiều hơn về điều đó sớm).

f = 10  # Frequency, in cycles per second, or Hertz
f_s = 100  # Sampling rate, or number of measurements per second

t = np.linspace(0, 2, 2 * f_s, endpoint=False)
x = np.sin(f * 2 * np.pi * t)

fig, ax = plt.subplots()
ax.plot(t, x)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Signal amplitude');

8

Trả lại tần số mẫu DFT.

f = 10  # Frequency, in cycles per second, or Hertz
f_s = 100  # Sampling rate, or number of measurements per second

t = np.linspace(0, 2, 2 * f_s, endpoint=False)
x = np.sin(f * 2 * np.pi * t)

fig, ax = plt.subplots()
ax.plot(t, x)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Signal amplitude');

9

Tính toán DFT của một chuỗi thực, khai thác tính đối xứng của phổ kết quả để tăng hiệu suất. Cũng được sử dụng bởi

import numpy as np

2 trong nội bộ khi áp dụng.

Danh sách này được bổ sung bởi các chức năng sau trong Numpy:

from scipy import fftpack

X = fftpack.fft(x)
freqs = fftpack.fftfreq(len(x)) * f_s

fig, ax = plt.subplots()

ax.stem(freqs, np.abs(X))
ax.set_xlabel('Frequency in Hertz [Hz]')
ax.set_ylabel('Frequency Domain (Spectrum) Magnitude')
ax.set_xlim(-f_s / 2, f_s / 2)
ax.set_ylim(-5, 110)

1,

from scipy import fftpack

X = fftpack.fft(x)
freqs = fftpack.fftfreq(len(x)) * f_s

fig, ax = plt.subplots()

ax.stem(freqs, np.abs(X))
ax.set_xlabel('Frequency in Hertz [Hz]')
ax.set_ylabel('Frequency Domain (Spectrum) Magnitude')
ax.set_xlim(-f_s / 2, f_s / 2)
ax.set_ylim(-5, 110)

2,

from scipy import fftpack

X = fftpack.fft(x)
freqs = fftpack.fftfreq(len(x)) * f_s

fig, ax = plt.subplots()

ax.stem(freqs, np.abs(X))
ax.set_xlabel('Frequency in Hertz [Hz]')
ax.set_ylabel('Frequency Domain (Spectrum) Magnitude')
ax.set_xlim(-f_s / 2, f_s / 2)
ax.set_ylim(-5, 110)

3,

from scipy import fftpack

X = fftpack.fft(x)
freqs = fftpack.fftfreq(len(x)) * f_s

fig, ax = plt.subplots()

ax.stem(freqs, np.abs(X))
ax.set_xlabel('Frequency in Hertz [Hz]')
ax.set_ylabel('Frequency Domain (Spectrum) Magnitude')
ax.set_xlim(-f_s / 2, f_s / 2)
ax.set_ylim(-5, 110)

4,

from scipy import fftpack

X = fftpack.fft(x)
freqs = fftpack.fftfreq(len(x)) * f_s

fig, ax = plt.subplots()

ax.stem(freqs, np.abs(X))
ax.set_xlabel('Frequency in Hertz [Hz]')
ax.set_ylabel('Frequency Domain (Spectrum) Magnitude')
ax.set_xlim(-f_s / 2, f_s / 2)
ax.set_ylim(-5, 110)

5

Chức năng tiêu gió thon.

DFT cũng được sử dụng để thực hiện các kết hợp nhanh của các đầu vào lớn bằng

from scipy import fftpack

X = fftpack.fft(x)
freqs = fftpack.fftfreq(len(x)) * f_s

fig, ax = plt.subplots()

ax.stem(freqs, np.abs(X))
ax.set_xlabel('Frequency in Hertz [Hz]')
ax.set_ylabel('Frequency Domain (Spectrum) Magnitude')
ax.set_xlim(-f_s / 2, f_s / 2)
ax.set_ylim(-5, 110)

6.

Scipy kết thúc thư viện FFTPack Fortran, đây không phải là loại nhanh nhất ngoài kia, nhưng không giống như các gói như FFTW, nó có giấy phép phần mềm miễn phí cho phép.

Chọn độ dài của DFT

Một tính toán ngây thơ của các hoạt động DFT (N2 ).6 Làm thế nào? Vâng, bạn có các hình sin n (phức tạp) có tần số khác nhau (2πf × 0, 2πf × 1; 2πf × 3, ..., 2πf × (n - 1)) và bạn muốn thấy tín hiệu của bạn tương ứng mạnh như thế nào với mỗi . Bắt đầu với lần đầu tiên, bạn lấy sản phẩm chấm với tín hiệu (bản thân nó, đòi hỏi n hoạt động nhân). Lặp lại hoạt động này n lần, một lần cho mỗi hình sin, sau đó đưa ra các hoạt động N2.

(N2) operations.6 How come? Well, you have N (complex) sinusoids of different frequencies (2πf × 0, 2πf × 1; 2πf × 3, ..., 2πf × (N – 1)), and you want to see how strongly your signal corresponds to each. Starting with the first, you take the dot-product with the signal (which, in itself, entails N multiplication operations). Repeating this operation N times, once for each sinusoid, then gives N2 operations.

Bây giờ, tương phản rằng với FFT, đó là (n log n) trong trường hợp lý tưởng do tái sử dụng thông minh các tính toán, một sự cải thiện lớn! Tuy nhiên, thuật toán Cooley-Tukey cổ điển được triển khai trong FFTPack (và được sử dụng bởi SCIPY) đã chia lại một cách biến đổi thành các mảnh nhỏ hơn (có kích thước nguyên tố) và chỉ hiển thị sự cải thiện này cho độ dài đầu vào mượt mà (một chiều dài đầu vào được coi là trơn tru khi nó Hệ số nguyên tố lớn nhất là nhỏ, như trong Hình & NBSP; 4-5). Đối với các mảnh có kích thước nguyên tố lớn, các thuật toán Bluestein hoặc Rader có thể được sử dụng cùng với thuật toán Cooley-Tukey, nhưng tối ưu hóa này không được thực hiện trong FFTPack.7

(N log N) in the ideal case due to the clever reuse of calculations—a great improvement! However, the classical Cooley-Tukey algorithm implemented in FFTPACK (and used by SciPy) recursively breaks up the transform into smaller (prime-sized) pieces and only shows this improvement for “smooth” input lengths (an input length is considered smooth when its largest prime factor is small, as shown in Figure 4-5). For large prime-sized pieces, the Bluestein or Rader algorithms can be used in conjunction with the Cooley-Tukey algorithm, but this optimization is not implemented in FFTPACK.7

Hãy để chúng tôi minh họa:

freq = fft.fftfreq(len(spectrum))
plot(freq, abs(spectrum))

3

Hình 4-5. Thời gian thực hiện FFT so với độ mịn cho các độ dài đầu vào khác nhauFFT execution time versus smoothness for different input lengths

Trực giác là, đối với số lượng mượt mà, FFT có thể được chia thành nhiều mảnh nhỏ. Sau khi thực hiện FFT trên phần đầu tiên, chúng ta có thể sử dụng lại các kết quả đó trong các tính toán tiếp theo. Điều này giải thích lý do tại sao chúng tôi đã chọn chiều dài 1.024 cho các lát âm thanh của chúng tôi trước đó, nó có độ mượt mà chỉ 2, dẫn đến thuật toán Cooley-2-Tukey tối ưu tối ưu, tính toán FFT chỉ bằng cách sử dụng (n/2) log2N = log2N = 5.120 phép nhân phức, thay vì N2 = 1.048.576. Chọn N = 2M luôn đảm bảo N mịn tối đa (và do đó, FFT nhanh nhất).

Nhiều khái niệm DFT hơn

Tiếp theo, chúng tôi trình bày một vài khái niệm phổ biến đáng để biết trước khi vận hành máy móc biến đổi Fourier nặng, sau đó chúng tôi giải quyết một vấn đề trong thế giới thực khác: phân tích phát hiện mục tiêu trong dữ liệu radar.

Tần số và thứ tự của chúng

Vì lý do lịch sử, hầu hết các triển khai trả lại một mảng trong đó các tần số thay đổi từ thấp đến cao đến thấp (xem Fourier rời rạc biến đổi để giải thích thêm về tần số). Ví dụ: khi chúng ta thực hiện biến đổi Fourier thực của tín hiệu của tất cả các tín hiệu, một đầu vào không có biến thể và do đó chỉ có thành phần Fourier không đổi chậm nhất Jargon cho giá trị trung bình của tín hiệu), xuất hiện dưới dạng mục đầu tiên:

freq = fft.fftfreq(len(spectrum))
plot(freq, abs(spectrum))

4

freq = fft.fftfreq(len(spectrum))
plot(freq, abs(spectrum))

5

Khi chúng tôi thử FFT trên tín hiệu thay đổi nhanh chóng, chúng tôi sẽ thấy một thành phần tần số cao xuất hiện:

freq = fft.fftfreq(len(spectrum))
plot(freq, abs(spectrum))

6

freq = fft.fftfreq(len(spectrum))
plot(freq, abs(spectrum))

7

freq = fft.fftfreq(len(spectrum))
plot(freq, abs(spectrum))

8

Lưu ý rằng FFT trả về một phổ phức tạp, trong trường hợp đầu vào thực, là đối xứng liên hợp (nghĩa là, đối xứng trong phần thực và đối xứng trong phần tưởng tượng):

freq = fft.fftfreq(len(spectrum))
plot(freq, abs(spectrum))

9

threshold = 0.5 * max(abs(spectrum))
mask = abs(spectrum) > threshold
peaks = freq[mask]

0

(Và, một lần nữa, nhớ lại rằng thành phần đầu tiên là

from scipy import fftpack

X = fftpack.fft(x)
freqs = fftpack.fftfreq(len(x)) * f_s

fig, ax = plt.subplots()

ax.stem(freqs, np.abs(X))
ax.set_xlabel('Frequency in Hertz [Hz]')
ax.set_ylabel('Frequency Domain (Spectrum) Magnitude')
ax.set_xlim(-f_s / 2, f_s / 2)
ax.set_ylim(-5, 110)

7.)

Hàm

f = 10  # Frequency, in cycles per second, or Hertz
f_s = 100  # Sampling rate, or number of measurements per second

t = np.linspace(0, 2, 2 * f_s, endpoint=False)
x = np.sin(f * 2 * np.pi * t)

fig, ax = plt.subplots()
ax.plot(t, x)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Signal amplitude');

8 cho chúng ta biết tần số nào chúng ta đang xem xét cụ thể:

threshold = 0.5 * max(abs(spectrum))
mask = abs(spectrum) > threshold
peaks = freq[mask]

1

threshold = 0.5 * max(abs(spectrum))
mask = abs(spectrum) > threshold
peaks = freq[mask]

2

Kết quả cho chúng ta biết rằng thành phần tối đa của chúng ta xảy ra ở tần số 0,5 chu kỳ trên mỗi mẫu. Điều này đồng ý với đầu vào, trong đó một chu kỳ trừ một lần lặp lại mỗi mẫu thứ hai.

Đôi khi, thật thuận tiện khi xem phổ được tổ chức hơi khác nhau, từ âm tính cao đến thấp đến thấp (hiện tại, chúng ta sẽ không đi sâu vào khái niệm về tần số âm, ngoài việc nói một Sóng sin thế giới được tạo ra bởi sự kết hợp của tần số dương và âm). Chúng tôi cải tổ phổ bằng hàm

f = 10  # Frequency, in cycles per second, or Hertz
f_s = 100  # Sampling rate, or number of measurements per second

t = np.linspace(0, 2, 2 * f_s, endpoint=False)
x = np.sin(f * 2 * np.pi * t)

fig, ax = plt.subplots()
ax.plot(t, x)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Signal amplitude');

6.

Hãy để kiểm tra các thành phần tần số trong hình ảnh ồn ào (Hình & NBSP; 4-7). Lưu ý rằng, trong khi một hình ảnh tĩnh không có thành phần thay đổi theo thời gian, các giá trị của nó thay đổi trên không gian. DFT áp dụng như nhau cho một trong hai trường hợp.

Đầu tiên, tải và hiển thị hình ảnh:

threshold = 0.5 * max(abs(spectrum))
mask = abs(spectrum) > threshold
peaks = freq[mask]

3

threshold = 0.5 * max(abs(spectrum))
mask = abs(spectrum) > threshold
peaks = freq[mask]

4

Hình 4-7. Một hình ảnh ồn ào của cuộc đổ bộ mặt trăngA noisy image of the moon landing

Không điều chỉnh màn hình của bạn! Hình ảnh bạn đang nhìn thấy là có thật, mặc dù rõ ràng bị biến dạng bởi các thiết bị đo hoặc truyền.

Để kiểm tra phổ của hình ảnh, chúng tôi sử dụng

import numpy as np

7 (thay vì

import numpy as np

2) để tính toán DFT, vì nó có nhiều hơn một chiều. FFT 2D tương đương với việc lấy FFT 1D trên các hàng và sau đó trên các cột hoặc ngược lại.

threshold = 0.5 * max(abs(spectrum))
mask = abs(spectrum) > threshold
peaks = freq[mask]

5

Một lần nữa, chúng tôi lấy nhật ký của phổ để nén phạm vi của các giá trị, trước khi hiển thị:

threshold = 0.5 * max(abs(spectrum))
mask = abs(spectrum) > threshold
peaks = freq[mask]

6

Lưu ý các giá trị cao xung quanh nguồn gốc (giữa) của quang phổ Các hệ số này mô tả các tần số thấp hoặc các phần mịn của hình ảnh, một bức tranh mơ hồ của ảnh. Các thành phần tần số cao hơn, trải khắp phổ, điền vào các cạnh và chi tiết. Các đỉnh xung quanh tần số cao hơn tương ứng với nhiễu định kỳ.

Từ ảnh, chúng ta có thể thấy rằng tiếng ồn (tạo tác đo lường) có tính định kỳ cao, vì vậy chúng tôi hy vọng sẽ loại bỏ nó bằng cách loại bỏ các phần tương ứng của phổ (Hình & NBSP; 4-8).

Hình ảnh với những đỉnh núi bị đàn áp thực sự trông khá khác biệt!

threshold = 0.5 * max(abs(spectrum))
mask = abs(spectrum) > threshold
peaks = freq[mask]

7

Hình 4-8. Hình ảnh hạ cánh mặt trăng được lọc và quang phổ của nóFiltered moon landing image and its spectrum

Cửa sổ

Nếu chúng ta kiểm tra biến đổi Fourier của xung hình chữ nhật, chúng ta sẽ thấy các thùy bên đáng kể trong quang phổ:

threshold = 0.5 * max(abs(spectrum))
mask = abs(spectrum) > threshold
peaks = freq[mask]

8

Về lý thuyết, bạn sẽ cần một sự kết hợp của nhiều hình sin (tần số) vô hạn để thể hiện bất kỳ sự chuyển đổi đột ngột nào; Các hệ số thường sẽ có cấu trúc thùy bên giống như được thấy ở đây cho xung.

Điều quan trọng, DFT giả định rằng tín hiệu đầu vào là định kỳ. Nếu tín hiệu không, giả định chỉ đơn giản là, ngay ở cuối tín hiệu, nó sẽ quay trở lại giá trị bắt đầu của nó. Xem xét hàm, x (t), hiển thị ở đây:

Chúng tôi đo tín hiệu chỉ trong một thời gian ngắn, được dán nhãn TEFF. Biến đổi Fourier giả định rằng x (8) = x (0) và tín hiệu được tiếp tục là đường nét, thay vì đường rắn,. Điều này giới thiệu một bước nhảy lớn ở rìa, với sự dao động dự kiến ​​trong quang phổ:

threshold = 0.5 * max(abs(spectrum))
mask = abs(spectrum) > threshold
peaks = freq[mask]

9

Thay vì hai dòng dự kiến, các đỉnh được trải ra trong quang phổ.

Chúng ta có thể chống lại hiệu ứng này bằng một quá trình gọi là cửa sổ. Hàm ban đầu được nhân với hàm cửa sổ như cửa sổ Kaiser K (n,). Ở đây chúng tôi trực quan hóa nó trong khoảng từ 0 đến 100:

t = arange(4e6) / 1e6 # sampling times in seconds
data = sin(2 * pi * 1000 * t)

0

Bằng cách thay đổi tham số, chúng ta có thể thay đổi hình dạng của cửa sổ từ hình chữ nhật (= 0, không có cửa sổ) thành một cửa sổ tạo ra các tín hiệu tăng trơn tru từ 0 và giảm xuống 0 tại các điểm cuối của khoảng thời gian được lấy mẫu, tạo ra rất thấp thùy bên (thường từ 5 đến 10) .9

Áp dụng cửa sổ Kaiser ở đây, chúng tôi thấy rằng các thùy bên đã giảm đáng kể, với chi phí mở rộng một chút trong thùy chính.

Hiệu quả của việc bán hàng ví dụ trước đây của chúng tôi là đáng chú ý:

t = arange(4e6) / 1e6 # sampling times in seconds
data = sin(2 * pi * 1000 * t)

1

Ứng dụng trong thế giới thực: Phân tích dữ liệu radar

Các radar FMCW được điều chế tuyến tính (sóng liên tục điều chế tần số) sử dụng rộng rãi thuật toán FFT để xử lý tín hiệu và cung cấp các ví dụ về các ứng dụng khác nhau của FFT. Chúng tôi sẽ sử dụng dữ liệu thực tế từ radar FMCW để chứng minh một ứng dụng như vậy: phát hiện mục tiêu.

Một cách thô bạo, một radar FMCW hoạt động như thế này (xem một hệ thống radar FMCW đơn giản và Hình & NBSP; 4-9 để biết thêm chi tiết):

1. Một tín hiệu với tần số thay đổi được tạo ra. Tín hiệu này được truyền bởi một ăng -ten, sau đó nó đi ra ngoài, ra khỏi radar. Khi nó chạm vào một vật thể, một phần của tín hiệu được phản xạ lại vào radar, nơi nó được nhận, nhân với một bản sao của tín hiệu truyền và lấy mẫu, biến nó thành các số được đóng gói thành một mảng. Thách thức của chúng tôi là giải thích những con số đó để hình thành kết quả có ý nghĩa.

2. Bước nhân trước là quan trọng. Từ trường học, hãy nhớ lại danh tính lượng giác:

3. Do đó, nếu chúng ta nhân tín hiệu nhận được bằng tín hiệu truyền, chúng ta hy vọng hai thành phần tần số sẽ xuất hiện trong phổ: một thành phần khác biệt về tần số giữa tín hiệu nhận và truyền và một phần là tổng của tần số của chúng.

4. Chúng tôi đặc biệt quan tâm đến lần đầu tiên, vì điều đó cho chúng tôi một số dấu hiệu cho thấy tín hiệu mất bao lâu để phản ánh trở lại radar (nói cách khác, đối tượng cách chúng tôi bao xa!). Chúng tôi loại bỏ cái khác bằng cách áp dụng bộ lọc thông thấp vào tín hiệu (nghĩa là, bộ lọc loại bỏ mọi tần số cao).

Tóm lại, chúng ta nên lưu ý rằng:

• Dữ liệu đến máy tính bao gồm N mẫu được lấy mẫu (từ tín hiệu được nhân, được lọc) ở tần số mẫu của FS.
• Biên độ của tín hiệu được trả về thay đổi tùy thuộc vào cường độ của phản xạ (nghĩa là, là một thuộc tính của đối tượng mục tiêu và khoảng cách giữa mục tiêu và radar).
• Tần số đo được là một dấu hiệu của khoảng cách của đối tượng mục tiêu từ radar.

Để bắt đầu phân tích dữ liệu radar, chúng tôi sẽ tạo ra một số tín hiệu tổng hợp, sau đó chúng tôi sẽ chuyển trọng tâm sang đầu ra của một radar thực tế.

Hãy nhớ lại rằng radar đang tăng tần số của nó khi nó truyền với tốc độ S Hz/s. Sau một khoảng thời gian nhất định, T, đã qua, tần số bây giờ sẽ cao hơn TS (Hình & NBSP; 4-10). Trong cùng khoảng thời gian đó, tín hiệu radar đã di chuyển d = t/v mét, trong đó v là tốc độ của sóng truyền qua không khí (gần giống với tốc độ của ánh sáng, 3 × 108 m/s).

Hình 4-10. Các mối quan hệ tần số trong radar FMCW với điều chế tần số tuyến tínhThe frequency relationships in an FMCW radar with linear frequency modulation

Kết hợp các quan sát trên, chúng ta có thể tính toán lượng thời gian sẽ đưa tín hiệu đi đến, bật ra và trở lại từ một mục tiêu cách xa R xa:

tr = 2r/v

t = arange(4e6) / 1e6 # sampling times in seconds
data = sin(2 * pi * 1000 * t)

2

Ở đây, chúng tôi đã tạo ra một tín hiệu tổng hợp, VSingle, nhận được khi nhìn vào một mục tiêu duy nhất (xem Hình & NBSP; 4-11). Bằng cách đếm số lượng chu kỳ nhìn thấy trong một khoảng thời gian nhất định, chúng ta có thể tính toán tần số của tín hiệu và do đó khoảng cách đến mục tiêu.

Một radar thực sự sẽ hiếm khi chỉ nhận được một tiếng vang duy nhất. Tín hiệu mô phỏng VSIM cho thấy tín hiệu radar sẽ trông như thế nào với năm mục tiêu ở các phạm vi khác nhau (bao gồm hai phạm vi gần nhau ở mức 154 và 159 mét) và Vactual (T) cho thấy tín hiệu đầu ra thu được bằng radar thực tế. Khi chúng ta thêm nhiều tiếng vang lại với nhau, kết quả có ý nghĩa trực quan ít (Hình & NBSP; 4-11); Cho đến khi, đó là, chúng tôi nhìn nó cẩn thận hơn qua lăng kính của DFT.

Hình 4-11. Tín hiệu đầu ra của máy thu: (a) mục tiêu mô phỏng đơn, (b) năm mục tiêu mô phỏng và (c) dữ liệu radar thực tếReceiver output signals: (a) single simulated target, (b) five simulated targets, and (c) actual radar data

Dữ liệu radar trong thế giới thực được đọc từ tệp .npz định dạng không định dạng (một định dạng lưu trữ tương thích chéo, đa nền tảng và đa nền tảng). Các tệp này có thể được lưu với các hàm

(-5, 110)

2 hoặc

(-5, 110)

3. Lưu ý rằng mô hình con SCIPY từ

(-5, 110)

4 cũng có thể dễ dàng đọc các định dạng khác, chẳng hạn như các tệp MATLAB và NETCDF.

t = arange(4e6) / 1e6 # sampling times in seconds
data = sin(2 * pi * 1000 * t)

3

Vì các tệp .npz có thể lưu trữ nhiều biến, chúng tôi phải chọn bản chúng tôi muốn:

(-5, 110)

5. Đó trả về một mảng numpy có cấu trúc với các trường sau:

(-5, 110)

6

Số nguyên 64 bit (8 byte) không dấu (

(-5, 110)

7)

(-5, 110)

8

Số nguyên không dấu 32-bit (4 byte) (

(-5, 110)

9)

import numpy.fft as fft
spectrum = fft.fft(data)

00

import numpy.fft as fft
spectrum = fft.fft(data)

01

Phao 32 bit (

import numpy.fft as fft
spectrum = fft.fft(data)

02)

import numpy.fft as fft
spectrum = fft.fft(data)

03

Phao 32 bit (

import numpy.fft as fft
spectrum = fft.fft(data)

02)

import numpy.fft as fft
spectrum = fft.fft(data)

05

Số nguyên 8 bit (1 byte) không dấu (

import numpy.fft as fft
spectrum = fft.fft(data)

06)

import numpy.fft as fft
spectrum = fft.fft(data)

07

Số nguyên 16 bit (2 byte) không dấu (

import numpy.fft as fft
spectrum = fft.fft(data)

08)

import numpy.fft as fft
spectrum = fft.fft(data)

09

Số nguyên 8 bit (1 byte) không dấu (

import numpy.fft as fft
spectrum = fft.fft(data)

10)

import numpy.fft as fft
spectrum = fft.fft(data)

11

2.048 số nguyên 16 bit (2 byte) không dấu (

import numpy.fft as fft
spectrum = fft.fft(data)

08)

Mặc dù đúng là các mảng numpy là đồng nhất (tức là, tất cả các yếu tố bên trong đều giống nhau), điều đó không có nghĩa là các yếu tố đó không thể là các yếu tố hợp chất, như trường hợp ở đây.

Một trường riêng lẻ được truy cập bằng cú pháp từ điển:

t = arange(4e6) / 1e6 # sampling times in seconds
data = sin(2 * pi * 1000 * t)

4

Để tóm tắt những gì chúng tôi đã thấy cho đến nay: các phép đo được hiển thị (VSIM và Vactual) là tổng các tín hiệu hình sin được phản ánh bởi mỗi đối tượng. Chúng ta cần xác định từng thành phần cấu thành của các tín hiệu radar tổng hợp này. FFT là công cụ sẽ làm điều này cho chúng tôi.

Thuộc tính tín hiệu trong miền tần số

Đầu tiên, chúng tôi lấy FFT của ba tín hiệu của chúng tôi (mục tiêu đơn tổng hợp, đa mục tiêu tổng hợp và thực) và sau đó hiển thị các thành phần tần số dương (nghĩa là, các thành phần 0 đến n/2; xem Hình & NBSP; 4-12). Chúng được gọi là dấu vết phạm vi trong thuật ngữ radar.

t = arange(4e6) / 1e6 # sampling times in seconds
data = sin(2 * pi * 1000 * t)

5

Hình 4-12. Dấu vết phạm vi cho: (a) mục tiêu mô phỏng đơn, (b) nhiều mục tiêu mô phỏng và (c) các mục tiêu trong thế giới thựcRange traces for: (a) single simulated target, (b) multiple simulated targets, and (c) real-world targets

Đột nhiên, thông tin có ý nghĩa!

Cốt truyện cho | v0 | hiển thị rõ ràng mục tiêu tại Thành phần 67 và cho | VSIM | hiển thị các mục tiêu tạo ra tín hiệu không thể giải thích được trong miền thời gian. Tín hiệu radar thực, | Vactual |, cho thấy một số lượng lớn mục tiêu giữa thành phần 400 và 500 với một đỉnh lớn trong thành phần 443. Đây là một sự trở lại của tiếng vang từ một radar chiếu sáng bức tường cao của mỏ mở.

Để có được thông tin hữu ích từ cốt truyện, chúng tôi phải xác định phạm vi! Một lần nữa, chúng tôi sử dụng công thức:

Trong thuật ngữ radar, mỗi thành phần DFT được gọi là thùng phạm vi.

Phương trình này cũng xác định độ phân giải phạm vi của radar: các mục tiêu sẽ chỉ có thể phân biệt được nếu chúng được phân tách bằng nhiều hơn hai thùng phạm vi, ví dụ::

Đây là một tài sản cơ bản của tất cả các loại radar.

Kết quả này khá thỏa mãn, nhưng phạm vi động lớn đến mức chúng ta có thể rất dễ dàng bỏ lỡ một số đỉnh. Hãy cùng lấy nhật ký như trước đây với phổ phổ:

t = arange(4e6) / 1e6 # sampling times in seconds
data = sin(2 * pi * 1000 * t)

6

Phạm vi động có thể quan sát được được cải thiện nhiều trong các ô này. Ví dụ, trong tín hiệu radar thực, sàn nhiễu của radar đã được nhìn thấy (nghĩa là, mức độ nhiễu điện tử trong hệ thống bắt đầu hạn chế khả năng phát hiện ra mục tiêu radar radar).

Cửa sổ, áp dụng

Chúng tôi đang đến đó, nhưng trong quang phổ của tín hiệu mô phỏng, chúng tôi vẫn không thể phân biệt các đỉnh ở mức 154 và 159 mét. Ai biết những gì chúng tôi thiếu trong tín hiệu trong thế giới thực! Để làm sắc nét các đỉnh, chúng tôi sẽ trở lại hộp công cụ của mình và sử dụng cửa sổ.

Dưới đây là các tín hiệu được sử dụng cho đến nay trong ví dụ này, được cửa sổ với cửa sổ Kaiser với β = 6.1:

t = arange(4e6) / 1e6 # sampling times in seconds
data = sin(2 * pi * 1000 * t)

7

Và các dấu vết FFTS, hoặc phạm vi tương ứng, theo thuật ngữ radar:

t = arange(4e6) / 1e6 # sampling times in seconds
data = sin(2 * pi * 1000 * t)

8

So sánh những điều này với các dấu vết phạm vi trước đó. Có một mức độ hạ thấp đáng kể về mức độ thù hận, nhưng ở mức giá: các đỉnh đã thay đổi hình dạng, mở rộng và trở nên ít đỉnh hơn, do đó hạ thấp độ phân giải radar, đó là khả năng của radar để phân biệt giữa hai mục tiêu cách nhau gần nhau . Sự lựa chọn của cửa sổ là một sự thỏa hiệp giữa mức độ và độ phân giải bên cạnh. Mặc dù vậy, đề cập đến dấu vết cho VSIM, cửa sổ đã tăng đáng kể khả năng phân biệt mục tiêu nhỏ với hàng xóm lớn của nó.

Trong dấu vết phạm vi dữ liệu radar thực, cửa sổ cũng đã làm giảm các thùy bên. Điều này có thể nhìn thấy nhất ở độ sâu của notch giữa hai nhóm mục tiêu.

Hình ảnh radar

Biết cách phân tích một dấu vết duy nhất, chúng ta có thể mở rộng để xem xét hình ảnh radar.

Dữ liệu được tạo ra bởi một radar với ăng -ten phản xạ parabol. Nó tạo ra một chùm bút chì tròn chỉ thị với góc lan rộng hai độ giữa các điểm nửa công suất. Khi được định hướng với tỷ lệ mắc bình thường tại một mặt phẳng, radar sẽ chiếu sáng một điểm có đường kính khoảng 2m ở khoảng cách 60m. Bên ngoài vị trí này, sức mạnh giảm khá nhanh, nhưng tiếng vang mạnh mẽ từ bên ngoài vị trí vẫn sẽ vẫn được nhìn thấy.

Bằng cách thay đổi Azimuth chùm bút chì (vị trí bên trái bên phải) và độ cao (vị trí lên xuống), chúng ta có thể quét nó qua khu vực mục tiêu quan tâm. Khi các phản xạ được chọn, chúng ta có thể tính khoảng cách đến gương phản xạ (đối tượng được đánh bởi tín hiệu radar). Cùng với vị thế và độ cao của chùm bút chì hiện tại, điều này xác định vị trí phản xạ trong 3D.

Một độ dốc đá bao gồm hàng ngàn phản xạ. Một thùng phạm vi có thể được coi là một quả cầu lớn với radar ở trung tâm của nó giao nhau với độ dốc dọc theo một đường rách. Các bộ phân tán trên dòng này sẽ tạo ra các phản xạ trong phạm vi này. Bước sóng của radar (khoảng cách sóng truyền di chuyển trong một lần dao động thứ hai) là khoảng 30 mm. Các phản xạ từ các bộ phân tán được phân tách bằng bội số lẻ của một bước sóng trong phạm vi, khoảng 7,5mm, sẽ có xu hướng can thiệp phá hủy, trong khi các chất phân tán được phân tách bởi bội số của một nửa bước sóng sẽ có xu hướng can thiệp xây dựng ở radar. Các phản xạ kết hợp để tạo ra các điểm rõ ràng của các phản xạ mạnh. Radar cụ thể này di chuyển ăng -ten của nó để quét các vùng nhỏ bao gồm 20 độ góc phương vị và thùng độ cao 30 độ được quét theo các bước 0,5 độ.

Bây giờ chúng tôi sẽ vẽ một số sơ đồ đường viền của dữ liệu radar kết quả. Tham khảo Hình & NBSP; 4-13 để xem các lát cắt khác nhau được thực hiện như thế nào. Một lát cắt đầu tiên ở phạm vi cố định cho thấy sức mạnh của tiếng vang chống lại độ cao và phương vị. Hai lát khác ở độ cao cố định và góc phương vị, tương ứng, hiển thị độ dốc (xem Hình 4-13 và 4-14). Cấu trúc bước của bức tường cao trong mỏ Opencast có thể nhìn thấy trong mặt phẳng phương vị.

Hình 4-13. Sơ đồ hiển thị góc phương vị, độ cao và phạm vi qua khối lượng dữ liệuDiagram showing azimuth, elevation, and range slices through data volume

t = arange(4e6) / 1e6 # sampling times in seconds
data = sin(2 * pi * 1000 * t)

9

Hình 4-14. Các sơ đồ đường viền của các dấu vết phạm vi dọc theo các trục khác nhau (xem Hình & NBSP; 4-13)Contour plots of range traces along various axes (see Figure 4-13)

Trực quan hóa 3D

Chúng ta cũng có thể trực quan hóa âm lượng theo ba chiều (Hình & NBSP; 4-15).

Trước tiên chúng tôi tính toán ArgMax (chỉ số của giá trị tối đa) theo hướng phạm vi. Điều này sẽ đưa ra một dấu hiệu của phạm vi mà chùm tia radar chạm vào độ dốc đá. Mỗi chỉ số ArgMax được chuyển đổi thành tọa độ 3D (độ cao-azimuth):

# Make plots appear inline, set custom plotting style
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('style/elegant.mplstyle')

0

Lấy các tọa độ này, chúng tôi chiếu chúng lên mặt phẳng phương vị phương vị (bằng cách bỏ tọa độ phạm vi) và thực hiện một tesselation delaunay. Tesselation trả về một tập hợp các chỉ số vào các tọa độ của chúng tôi xác định các hình tam giác (hoặc đơn giản). Mặc dù các hình tam giác được nói một cách nghiêm ngặt được định nghĩa trên các tọa độ dự kiến, chúng tôi sử dụng các tọa độ ban đầu của chúng tôi để tái thiết, do đó thêm lại thành phần phạm vi:

# Make plots appear inline, set custom plotting style
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('style/elegant.mplstyle')

1

Đối với mục đích hiển thị, chúng tôi trao đổi trục phạm vi là đầu tiên:

# Make plots appear inline, set custom plotting style
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('style/elegant.mplstyle')

2

Bây giờ, Matplotlib từ

import numpy.fft as fft
spectrum = fft.fft(data)

13 có thể được sử dụng để trực quan hóa kết quả:

# Make plots appear inline, set custom plotting style
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('style/elegant.mplstyle')

3

Hình 4-15. Trực quan hóa 3D về vị trí độ dốc đá ước tính3D visualization of estimated rock slope position

Các ứng dụng tiếp theo của FFT

Các ví dụ trước chỉ cho thấy một trong những cách sử dụng của FFT trong radar. Có nhiều người khác, chẳng hạn như đo lường (Doppler) và nhận dạng mục tiêu. FFT có sức lan tỏa, và được nhìn thấy ở khắp mọi nơi từ MRI đến thống kê. Với các kỹ thuật cơ bản mà chương này phác thảo trong tay, bạn nên được trang bị tốt để sử dụng nó!

Đọc thêm

Về biến đổi Fourier:

• Athanasios Papoulis, Integral Fourier và các ứng dụng của nó (New York: McGraw-Hill, 1960).
• Ronald A. Bracewell, Biến đổi Fourier và các ứng dụng của nó (New York: McGraw-Hill, 1986).

Trên xử lý tín hiệu radar:

• Mark A. Richards, James A. Scheer và William A. Holm, Eds., Nguyên tắc của radar hiện đại: các nguyên tắc cơ bản (Raleigh, NC: Scitech, 2010).
• Mark A. Richards, Nguyên tắc cơ bản của xử lý tín hiệu radar (New York: McGraw-Hill, 2014).

Tập thể dục: Con chập hình ảnh

FFT thường được sử dụng để tăng tốc độ tích chập hình ảnh (chập là ứng dụng của bộ lọc trượt). Kết hợp một hình ảnh với

import numpy.fft as fft
spectrum = fft.fft(data)

14, sử dụng a) numpy từ

import numpy.fft as fft
spectrum = fft.fft(data)

15 và b)

import numpy.fft as fft
spectrum = fft.fft(data)

16. Xác nhận rằng kết quả là giống hệt nhau.

Hints:

• Việc tích chập

  import numpy.fft as fft
  spectrum = fft.fft(data)
  

  17 và

  import numpy.fft as fft
  spectrum = fft.fft(data)
  

  18 tương đương với

  import numpy.fft as fft
  spectrum = fft.fft(data)
  

  19, trong đó

  import numpy.fft as fft
  spectrum = fft.fft(data)
  

  20 và

  import numpy.fft as fft
  spectrum = fft.fft(data)
  

  21 lần lượt là FFT của X và Y.
• Để nhân

  import numpy.fft as fft
  spectrum = fft.fft(data)
  

  20 và

  import numpy.fft as fft
  spectrum = fft.fft(data)
  

  21, chúng phải có cùng kích thước. Sử dụng

  import numpy.fft as fft
  spectrum = fft.fft(data)
  

  24 để mở rộng

  import numpy.fft as fft
  spectrum = fft.fft(data)
  

  17 và

  import numpy.fft as fft
  spectrum = fft.fft(data)
  

  18 với số không (về phía bên phải và dưới) trước khi lấy FFT của họ.
• Bạn có thể thấy một số hiệu ứng cạnh. Bạn có thể loại bỏ chúng bằng cách tăng kích thước đệm, để cả

  import numpy.fft as fft
  spectrum = fft.fft(data)
  

  17 và

  import numpy.fft as fft
  spectrum = fft.fft(data)
  

  18 đều có kích thước

  import numpy.fft as fft
  spectrum = fft.fft(data)
  

  29.

Kiểm tra giải pháp trên mạng: hình ảnh chập hình ảnh.

1 DFT hoạt động trên dữ liệu được lấy mẫu, trái ngược với biến đổi Fourier tiêu chuẩn, được xác định cho các hàm liên tục.

2 Biến đổi Fourier về cơ bản cho chúng ta biết cách kết hợp một tập hợp các hình sin có tần số khác nhau để tạo thành tín hiệu đầu vào. Phổ bao gồm các số phức tạp, một cho mỗi hình sin. Một số phức tạp mã hóa hai điều: một độ lớn và một góc. Độ lớn là cường độ của xoang trong tín hiệu và góc là bao nhiêu nó được thay đổi trong thời gian. Tại thời điểm này, chúng tôi chỉ quan tâm đến độ lớn mà chúng tôi tính toán bằng cách sử dụng

import numpy.fft as fft
spectrum = fft.fft(data)

30.

3 Để biết thêm về các kỹ thuật để tính toán cả tần số (gần đúng) và thời gian xảy ra, hãy đọc lên phân tích sóng con.

4 Scipy đi đến một số nỗ lực để bảo tồn năng lượng trong quang phổ. Do đó, khi chỉ lấy một nửa các thành phần (với n chẵn), nó sẽ nhân các thành phần còn lại, ngoài các thành phần đầu tiên và cuối cùng, với hai thành phần đó (hai thành phần đó được chia sẻ bởi hai nửa của quang phổ). Nó cũng bình thường hóa cửa sổ bằng cách chia nó cho tổng của nó.

5 Thời gian trên thực tế, cũng có thể là vô hạn! Biến đổi Fourier liên tục chung cung cấp cho khả năng này. DFT thường được xác định trong một khoảng thời gian hữu hạn và đây là khoảng thời gian của hàm miền thời gian được chuyển đổi. Nói cách khác, nếu bạn thực hiện DFT nghịch đảo, bạn luôn nhận được tín hiệu định kỳ.

6 Trong khoa học máy tính, chi phí tính toán của một thuật toán thường được thể hiện trong ký hiệu Big Big O. Ký hiệu cho chúng ta một dấu hiệu cho thấy cách thức thực hiện thời gian thực hiện thuật toán với số lượng phần tử ngày càng tăng. Nếu một thuật toán là (n), thời gian thực hiện của nó tăng tuyến tính với số lượng các phần tử đầu vào (ví dụ: tìm kiếm một giá trị nhất định trong một danh sách chưa được phân loại là (n). Sắp xếp bong bóng là một ví dụ về thuật toán O (n2); Số lượng hoạt động chính xác được thực hiện, theo giả thuyết, có nghĩa là, có nghĩa là chi phí tính toán tăng theo bậc hai với số lượng các yếu tố đầu vào.

(N), its execution time increases linearly with the number of input elements (e.g., searching for a given value in an unsorted list is (N). Bubble sort is an example of an O(N2) algorithm; the exact number of operations performed may, hypothetically, be , meaning that the computational cost grows quadratically with the number of input elements.

7 Mặc dù lý tưởng nhất là chúng ta không muốn thực hiện lại các thuật toán hiện có, đôi khi nó trở nên cần thiết để có được tốc độ thực hiện tốt nhất có thể, và các công cụ như Cython, biên soạn Python với C và Numba, công cụ biên dịch đúng lúc Mã Python, làm cho cuộc sống dễ dàng hơn rất nhiều (và nhanh hơn!). Nếu bạn có thể sử dụng phần mềm được cấp phép GPL, bạn có thể xem xét sử dụng pyfftw cho các FFT nhanh hơn.

8 Chúng tôi để nó như một bài tập cho người đọc hình dung tình huống cho n lẻ. Trong chương này, tất cả các ví dụ sử dụng DFT đều theo thứ tự.

9 Các chức năng gió cổ điển bao gồm Hann, Hamming và Blackman. Chúng khác nhau về mức độ thùy bên của chúng và trong việc mở rộng thùy chính (trong miền Fourier). Một chức năng cửa sổ hiện đại và linh hoạt gần với tối ưu cho hầu hết các ứng dụng là cửa sổ Kaiser, một xấp xỉ tốt với cửa sổ hình cầu sinh học tối ưu, tập trung nhiều năng lượng nhất vào thùy chính. Chúng ta có thể điều chỉnh cửa sổ Kaiser để phù hợp với ứng dụng cụ thể, như được minh họa trong văn bản chính, bằng cách điều chỉnh tham số.

Làm thế nào để bạn tìm thấy tần số của tín hiệu trong Python?

Làm thế nào để bạn tìm thấy tần số của tín hiệu?

Làm thế nào để bạn vẽ một phổ tần số của tín hiệu trong Python?

FFTFREQ trở lại điều gì?