FM with IQ Demodulation
20161210[ Radio FM Math Demodulation ]
After some time reading about hacking RF signals:
 Using RTLSDR to automatically receive weather satellite images.
 Reverse engineering traffic lights with software defined radio.
 Spectrum painting on 2.4 GHz.
 RTLSDR Tutorial: Cheap ADSB aircraft radar
I bought a RTL2832U receiver, that can operate on a range of 701700 Mhz.
Libraries and Programs
After some time (6 months, hell yeah Brazil...) to receive the package, I begun to install all the programs and libraries to use the RTL.
$ sudo aptget install gnuradio rtlsdr gqrxsdr

gnuradio : GNU Radio is a free software development toolkit that provides signal processing blocks to implement softwaredefined radios and signal processing systems. It can be used with external RF hardware to create softwaredefined radios, or without hardware in a simulationlike environment. It is widely used in hobbyist, academic, and commercial environments to support both wireless communications research and realworld radio systems. (ty wikipedia).

rtlsdr have:
`
* rtl_adsb: a simple ADSB decoder for RTL2832 based DVBT receivers * rtl_eeprom: an EEPROM programming tool for RTL2832 based DVBT receivers * rtl_fm: a narrow band FM demodulator for RTL2832 based DVBT receivers * rtl_sdr: an I/Q recorder for RTL2832 based DVBT receivers * rtl_tcp: an I/Q spectrum server for RTL2832 based DVBT receivers * rtl_test: a benchmark tool for RTL2832 based DVBT receivers
`
 gqrxsdr is a RF receive interface that have a displays FFT plot and spectrum waterfall very helpful to detect signals, the software includes AM, SSB, FMN, FMW (mono and stereo) demodulators and Special FM mode for NOAA APT.
After installed all the software and plugging the receiver on the computer, we can see that have a Realtek chip with the lsusb feedback.
Bus 001 Device 004: ID 0bda:2838 Realtek Semiconductor Corp. RTL2838 DVBT
Ok, now with everything done, we can capture some data ! We can use the rtlsdr toolkit.
rtl_sdr capture.bin s 1.8e6 f 100.9e6
The capture.bin is the file that will have the IQ data, 1.8 MHz is the sample frequency and the frequency that we will look for is 100.9MHz (Rádio Atlântida), the radio station.
Warning !! 1.8e6 samples for second, so take 3~10 seconds because the fast growing of the file.
To demodulate we need do understand the modulation.
The FM is a method to modulate a message $m(t)$ with the change of the carrier ![cos(2\pi f_c t)] phase's . Wikipedia is your friend.
The phase have a relationship with the message, that is . So the simple FM modulator will look like:
The data will be transmitted after the modulation.
Demodulate FM, IQ Decomposition.
The RTL driver give to us the result of the QAM demodulator.
Where:
You can read more about IQ data here.
Now we know how to demodulate the FM, with the equations above, the result is:
Some code.
Now, with the knowledge of how the FM modulation and demodulation work, we can write a octave code that can show to us how the theory work in the real world.
%read file and convert data fid = fopen('./capture2.bin','rb'); y = fread(fid,'uint8=&amp;amp;amp;amp;amp;gt;double'); y = y127; % transform IQ data to an imaginary number yi = y(1:2:end)+i*y(2:2:end); % the radio frequency and the sample frequency freq = 100.9e6; fs = 1.8e6; % some math yang = angle(yi); % take the angle pi to pi yrap = unwrap(yang); % make the correction to continue the angle after pi and pi tdev = diff(yrap); % the derivative of theta sound(tdev,fs) % give us some music
You can download the capture.bin that I have used and the .m code to octave (open source rules) here: fm_demodulator_mono.
Stereo demodulation.
Ok, we have done the simpliest demodulation possible, without filters and a lot of things, so now we will use more energy to obtain a very good good
cool sound.
The FM spectrum
The FM spectrum have more things that I thought.
 30~15kHz L+R (mono), that we have already demodulated.
 19kHz, the omnipotent carrier.
 23~53kHz, a kind of dsbsc (double side band suppressedcarrier) with LR.
 55.35~58.65kHz, the RBDS (
RebeldesRadio Data system), a kind of digital data system.  58.65~76.65kHz, FTF first read this, now ask yourself wth Microsoft is here ??
 92+, IDKWII.
What is important ?
The FM modulation appears to be something like this:
Mathematically:
Where .
Yeah, this is a kind of crazy. The IQ data give to us the imaginary part that compose the , theoretically making the FFT of we can see something like the FM spectrum image.
% fft tdfft = fft(tdev); P2=abs(tdfft/ysize); P1=P2(1:ysize/2+1); P1(2:end1)=2*P1(2:end1); f=fs*(0:(ysize/2))/ysize; % plot the fft until 59kHz freqe=round(2*59e3*size(f)/fs); freqe=freqe(2); plot(f(1:freqe),(P1(1:freqe))) grid on xlabel 'freq Hz',ylabel 'dB' hold on axis([30 59e3 0 0.02]) Gab1 = [30 0; 30 0.02; 15e3 0.02; 15e3 0]; % R+L plot(Gab1(:,1),Gab1(:,2),'g'); Gab1 = [18.5e3 0; 18.5e3 0.02; 19.5e3 0.02; 19.5e3 0]; % carrier plot(Gab1(:,1),Gab1(:,2),'b'); Gab1 = [23e3 0; 23e3 0.02; 53e3 0.02; 53e3 0]; % RL plot(Gab1(:,1),Gab1(:,2),'r'); Gab1 = [55.35e3 0; 55.35e3 0.02; 58.65e3 0.02; 58.65e3 0]; % RBDS plot(Gab1(:,1),Gab1(:,2),'c');
And the result:
All good, all good sir !
The carrier and his friends.
Like said before, we have the IQ data that can give to us the and with a derivative we can have . Applying a lowpass filter on the L+R, a bandpass filter on the carrier, LR and RBDS we can have a clear signal of each data.
With the frequencies specifications, the filter was done:
%R+L [b,a] = butter(10,2*16e3/fs,'low'); R_p_L = filter(b,a,tdev); hhhh=R_p_L; plote %carrier [b,a] = butter(4,[2*18.5e3,2*19.5e3]/fs,'bandpass'); carrier = filter(b,a,tdev); hhhh=carrier; plote % RL [b,a] = butter(4,[2*23e3,2*53e3]/fs,'bandpass'); R_l_L = filter(b,a,tdev); hhhh=R_l_L; plote %RBDS [b,a] = butter(4,[2*55.35e3,2*58.65e3]/fs,'bandpass'); rbds = filter(b,a,tdev); hhhh=rbds; plote
The 'plote' code plot the fft of the signal.
To be continue..