DSP_experiments/FM_radio_sim.py at main · linux81/DSP_experiments · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert
from scipy.io.wavfile import write

# --- 1. Parametry symulacji ---
sampling_rate = 44100  # Hz
carrier_frequency = 100000  # Hz
left_audio_frequency = 440  # Hz (A4 note)
right_audio_frequency = 660  # Hz (E5 note)
pilot_frequency = 19000  # Hz
subcarrier_frequency = 38000  # Hz
modulation_index = 50 # Controls the frequency deviation

# --- 2. Generowanie sygnałów audio ---
t = np.linspace(0, 1, sampling_rate, endpoint=False)
left_signal = np.sin(2 * np.pi * left_audio_frequency * t)
right_signal = np.sin(2 * np.pi * right_audio_frequency * t)

# --- 3. Multipleksowanie sygnału stereo ---
mono_signal = left_signal + right_signal
difference_signal = left_signal - right_signal
subcarrier_signal = difference_signal * np.cos(2 * np.pi * subcarrier_frequency * t)
pilot_signal = np.cos(2 * np.pi * pilot_frequency * t)
multiplex_signal = mono_signal + subcarrier_signal + pilot_signal

# --- 4. Modulacja FM ---
instantaneous_phase = 2 * np.pi * carrier_frequency * t + modulation_index * np.cumsum(multiplex_signal) / sampling_rate
fm_signal = np.cos(instantaneous_phase)

# --- 5. Demodulacja FM (demodulator) ---
analytic_signal = hilbert(fm_signal)
instantaneous_phase = np.unwrap(np.angle(analytic_signal))
instantaneous_frequency = (np.diff(instantaneous_phase) / (2 * np.pi) * sampling_rate)
demodulated_multiplex_signal = instantaneous_frequency - carrier_frequency

# --- Normalizacja i usunięcie DC offset ---
demodulated_multiplex_signal -= np.mean(demodulated_multiplex_signal)
demodulated_multiplex_signal /= np.max(np.abs(demodulated_multiplex_signal))

# --- 6. Demultipleksowanie sygnału stereo ---
# Użycie odpowiedniej wielkości tablicy 't' dla demodulacji
t_demod = t[:-1]
demodulated_difference = demodulated_multiplex_signal * np.cos(2 * np.pi * subcarrier_frequency * t_demod)
demodulated_difference -= np.mean(demodulated_difference)
demodulated_difference /= np.max(np.abs(demodulated_difference))

# Odzyskiwanie sygnałów L i R z (L+R) i (L-R)
recovered_left = (demodulated_multiplex_signal + demodulated_difference) / 2
recovered_right = (demodulated_multiplex_signal - demodulated_difference) / 2

# --- 7. Zapis do pliku WAV ---
recovered_stereo = np.vstack((recovered_left, recovered_right)).T
recovered_stereo_int16 = np.int16(recovered_stereo * 32767)
write("demodulated_stereo_audio.wav", sampling_rate, recovered_stereo_int16)
print("Zdemodulowany sygnał stereo zapisano do 'demodulated_stereo_audio.wav'")

# --- 8. Wizualizacja przebiegów ---
plt.figure(figsize=(18, 12))
plt.subplot(3, 2, 1)
plt.title("Sygnały wejściowe (L i R)")
plt.plot(t[:1000], left_signal[:1000], label="Lewy (440 Hz)")
plt.plot(t[:1000], right_signal[:1000], label="Prawy (660 Hz)")
plt.xlabel("Czas (s)")
plt.ylabel("Amplituda")
plt.legend()
plt.grid(True)

plt.subplot(3, 2, 2)
plt.title("Złożony sygnał multipleksowy")
plt.plot(t[:1000], multiplex_signal[:1000])
plt.xlabel("Czas (s)")
plt.ylabel("Amplituda")
plt.grid(True)

plt.subplot(3, 2, 3)
plt.title("Sygnał FM (zmodulowany)")
plt.plot(t[:1000], fm_signal[:1000])
plt.xlabel("Czas (s)")
plt.ylabel("Amplituda")
plt.grid(True)

plt.subplot(3, 2, 4)
plt.title("Sygnał multipleksowy (po demodulacji)")
plt.plot(t_demod[:1000], demodulated_multiplex_signal[:1000])
plt.xlabel("Czas (s)")
plt.ylabel("Amplituda")
plt.grid(True)

plt.subplot(3, 2, 5)
plt.title("Odzyskany sygnał lewego kanału")
plt.plot(t_demod[:1000], recovered_left[:1000], label="Lewy (440 Hz)")
plt.xlabel("Czas (s)")
plt.ylabel("Amplituda")
plt.grid(True)

plt.subplot(3, 2, 6)
plt.title("Odzyskany sygnał prawego kanału")
plt.plot(t_demod[:1000], recovered_right[:1000], label="Prawy (660 Hz)")
plt.xlabel("Czas (s)")
plt.ylabel("Amplituda")
plt.grid(True)

plt.tight_layout()
plt.show()

# --- 9. Analiza FFT ---
def plot_fft(signal, title, fs, is_stereo=False):
    if is_stereo:
        signal = signal[:, 0]  # Take only the left channel for simplicity
    n = len(signal)
    fft_result = np.fft.fft(signal)
    fft_freq = np.fft.fftfreq(n, 1/fs)
    plt.plot(fft_freq[:n//2], np.abs(fft_result[:n//2]))
    plt.title(title)
    plt.xlabel("Częstotliwość (Hz)")
    plt.ylabel("Amplituda")
    plt.grid(True)

plt.figure(figsize=(18, 6))

plt.subplot(1, 3, 1)
plot_fft(multiplex_signal, "FFT sygnału multipleksowego", sampling_rate)

plt.subplot(1, 3, 2)
plot_fft(demodulated_multiplex_signal, "FFT po demodulacji", sampling_rate)
plt.xlim(0, 2500) # Ograniczenie widoku do sygnału audio

plt.subplot(1, 3, 3)
plot_fft(recovered_left, "FFT odzyskanego sygnału (L)", sampling_rate)
plt.xlim(0, 1000)

plt.tight_layout()
plt.show()