Pull to refresh

Кодируем и декодируем сообщение для внеземных цивилизаций

PythonProgrammingPopular scienceAstronauticsAstronomy
Привет Хабр.

Мотивом для этой статьи, на самом деле, послужил грустный повод. Всемирно известный радиотелескоп обсерватории Аресибо в Пуэрто-Рико разрушился и восстановлению не подлежит. Многие годы это был самый крупный радиотелескоп в мире с полной апертурой (диаметр 304 м, диапазон частот до 10 ГГц), с помощью которого было сделано немало открытий. Здесь на фото из Википедии он в еще рабочем состоянии:


Источник: en.wikipedia.org/wiki/Arecibo_Observatory

Но текст собственно о другом событии. В 1974 г с этого телескопа было отправлено в космос послание внеземным цивилизациям. Что и как было в нем закодировано, подробности под катом.

Кодирование


Для начала интересно понять, как сообщение передавалось. Как известно, размер сообщения был всего лишь 1679 бит (примерно 210 байт), и передавалось оно на частоте 2.38 ГГц с мощностью 450 КВт. Для передачи использовалась частотная модуляция со скоростью 10 бит/с. Число 1679 выбрано не случайно — это произведение двух простых чисел 23 и 73, таким образом, разложить картинку в виде прямоугольника можно единственным способом.

Готового радиосообщения в формате WAV мне найти не удалось, но в бинарном виде оно есть, а с помощью Python звук несложно сгенерировать самостоятельно. Желающие послушать, что услышат инопланетяне, могут скачать и запустить нижеприведенный код, который сгенерирует WAV файл. Для большей достоверности в сообщение также добавлен шум.

generate.py
import scipy.io.wavfile as wav
import scipy.signal as signal
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import os


message = """0000001010101000000000000101000001010000000100100010001000100
             1011001010101010101010100100100000000000000000000000000000000
             0000011000000000000000000011010000000000000000000110100000000
             0000000000101010000000000000000001111100000000000000000000000
             0000000001100001110001100001100010000000000000110010000110100
             0110001100001101011111011111011111011111000000000000000000000
             0000010000000000000000010000000000000000000000000000100000000
             0000000001111110000000000000111110000000000000000000000011000
             0110000111000110001000000010000000001000011010000110001110011
             0101111101111101111101111100000000000000000000000000100000011
             0000000001000000000001100000000000000010000011000000000011111
             1000001100000011111000000000011000000000000010000000010000000
             0100000100000011000000010000000110000110000001000000000011000
             1000011000000000000000110011000000000000011000100001100000000
             0110000110000001000000010000001000000001000001000000011000000
             0010001000000001100000000100010000000001000000010000010000000
             1000000010000000100000000000011000000000110000000011000000000
             1000111010110000000000010000000100000000000000100000111110000
             0000000010000101110100101101100000010011100100111111101110000
             1110000011011100000000010100000111011001000000101000001111110
             0100000010100000110000001000001101100000000000000000000000000
             0000000001110000010000000000000011101010001010101010100111000
             0000001010101000000000000000010100000000000000111110000000000
             0000001111111110000000000001110000000111000000000110000000000
             0110000000110100000000010110000011001100000001100110000100010
             1000001010001000010001001000100100010000000010001010001000000
             0000001000010000100000000000010000000001000000000000001001010
             00000000001111001111101001111000"""

def fftnoise(f):
    f = np.array(f, dtype='complex')
    n_p = (len(f) - 1) // 2
    phases = np.random.rand(n_p) * 2 * np.pi
    phases = np.cos(phases) + 1j * np.sin(phases)
    f[1:n_p+1] *= phases
    f[-1:-1-n_p:-1] = np.conj(f[1:n_p+1])
    return np.fft.ifft(f).real

def band_limited_noise(min_freq, max_freq, samples, samplerate=1):
    freqs = np.abs(np.fft.fftfreq(samples, 1/samplerate))
    f = np.zeros(samples)
    idx = np.where(np.logical_and(freqs>=min_freq, freqs<=max_freq))[0]
    f[idx] = 1
    return fftnoise(f)


message = ''.join(i for i in message if i.isdigit())
print("Original message:")
print(message)
print()

# Generate message
fs = 11025
f1, f2 = 3000, 4000
t_sym = 0.1
data = np.zeros(int(fs * t_sym * len(message)))
for p in range(len(message)):
    samples = np.linspace(0, t_sym, int(fs * t_sym), endpoint=False)
    freq = f2 if message[p] == '1' else f1
    data[int(fs * t_sym)*p:int(fs * t_sym)*(p + 1)] = 10000*(0.25*np.sin(2 * np.pi * freq * samples) + band_limited_noise(50, 5000, len(samples), fs))

wav.write('arecibo.wav', fs, np.int16(data))
print("WAV file saved")

Для удобства прослушивания я увеличил разнос частот, в оригинальном сообщении он был лишь 10 Гц. К сожалению, habr не позволяет приаттачивать звуковые файлы, желающие могут либо сгенерировать файл самостоятельно, либо использовать временную ссылку.

Да кстати, отправлено сообщение было в 1974 году. Вот туда:


Источник: en.wikipedia.org/wiki/Messier_13

Красивое звездное скопление M13 в созвездии Геркулеса, хорошо известное всем любителям астрономии, и доступное для наблюдений даже в небольшие телескопы. Скопление находится в 22 тыс световых лет, так что идти сообщение будет еще долго…

С кодированием мы разобрались, теперь представим, что мы приняли подобное сообщение — посмотрим как его можно декодировать.

Декодирование


Принцип частотной модуляции сам по себе прост — нулю и единице соответствуют разные частоты. На спектре это выглядит примерно так:



Существуют разные способы декодирования FSK, в качестве самого простого метода просто отфильтруем одну из частот:

fs, data = wav.read('arecibo.wav')

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = signal.butter(order, [low, high], btype='band')
    return b, a

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = signal.lfilter(b, a, data)
    return y

f1, f2 = 3000, 4000
data_f2 = butter_bandpass_filter(data, f2 - 200, f2 + 200, fs, order=3)

plt.plot(data)
plt.plot(data_f2)
plt.xlabel("Time")
plt.ylabel("Amplitude")
plt.title("Signal")
plt.show()

Результат нас вполне устраивает:



Разумеется, сигнал прошедший 22 тыс лет в космосе, скорее всего будет немного ослаблен, но для простоты будем считать что у инопланетян хорошие и не китайские радиоприемники…

Из картинки легко можно определить ширину одного бита. Нам нужно вывести биты в виде изображения. Т.к. послание отправлялось внеземной цивилизации — тем кто, по определению, не знает «земных» систем кодирования — передать растровую картинку было единственным логичным решением. В созвездии Геркулеса скорее всего не знают что такое ASCII или упаси бог, Unicode, но вывести растр на экран скорее всего, можно в любой точке Галактики. По крайней мере, цивилизация способная принять цифровой сигнал, скорее всего имеет и какое-то подобие монитора для его отображения.

Мы не знаем размера картинки, зато знаем размер одного бита и знаем размер всего сообщения. Можно просто перебрать все возможные варианты, благо их не так уж много:

ss = 1102  # Width of one symbol in samples
for iw in range(12*ss, 25*ss, ss):
    w, h = iw, 80
    image = Image.new('RGB', (w, h))

    px, py = 0, 0
    for p in range(data_f2.shape[0]):
        image.putpixel((px, py), (0, int(data_f2[p]//32), 0))
        px += 1
        if px >= w:
            px = 0
            py += 1
            if py >= h:
                break

    image = image.resize((w//10, 100*h))
    image.save("1/image-%d.png" % iw)


Для наглядности картинку пришлось растянуть, ибо 23 пиксела ширины по сегодняшним меркам это мягко говоря, маловато. Конечный результат вполне виден:



Окончательный вариант:



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

В картинке закодировано много всего (условно конечно), например вертикальная прямая шириной 2 пиксела над головой человека — это спираль ДНК (ведь очевидно же, разве нет?). Расшифровку остальных пиктограмм можно посмотреть по ссылке выше в Википедии.

Заключение


Как можно видеть, в 210 байт можно закодировать довольно много информации. Вообще, сама по себе задача отправки сигнала в дальний космос является далеко не простой, т.к. можно надеяться лишь на самые простые методы модуляции. Дойдет ли сообщение до адресата? Разумеется, скорее всего, вряд ли. Я не знаю, оценивали ли при отправке энергетику такой «линии связи» и примерную чувствительность приемника, нужную чтобы получить сигнал. Да это на самом деле, не так и важно — если подобные акции вдохновят кого-либо на изучение космоса, значит это было не зря. Ну а точный ответ мы сможем получить через 44 тыс лет, и я постараюсь обновить текст по мере поступления новых данных ;)
Tags:pythonSETIпосланиеНЛО
Hubs: Python Programming Popular science Astronautics Astronomy
Rating +26
Views 7k Add to bookmarks 31
Comments
Comments 35

Popular right now

Факультет Python-разработки
March 12, 2021180,000 ₽GeekBrains
Python QA Engineer
March 16, 202160,000 ₽OTUS
Python-разработчик с нуля
March 22, 202169,900 ₽Нетология
Python для работы с данными
March 26, 202131,500 ₽Нетология

Top of the last 24 hours