最近找到点时间继续开发sstv接收功能。 我想了一下，之前直接做martin m1这种彩色的sstv接收可能太复杂了，也许黑白的会简单一点，我在portapack上做的模拟电视解调也是黑白的，虽然图像不如彩色的好看，但是聊胜于无吧。所以我打算解调robot8BW。
现在接触下来fm解调的算法有4种，滤波器方法、FFT方法、希尔伯特变换、过零检测、相乘。
如果用第一种滤波器方法，只能解固定频率的信号，用来做vis解码还行，解图像肯定不行。
FFT从原来上好理解，但是如果不调用库，而是直接写代码实现比较复杂，而且需要一定长度的点来做运算，可能实时性不够。
希尔伯特变换和FFT有一样的问题，不调用现成的库自己从头写代码比较难。
过零检测是一种很原始的方法，早期sstv接收就是这样做的，通过数正弦波经过x轴的时间间隔来测量频率，理论和实现都简单，效果不太好。
相乘的方法自己实现也比较简单，效果也不错，毕竟之前已经在电脑上用它实现了martin m1的解调了，现在只是想把它用于黑白图片解调。
下面的一段代码，我用了过零检测、希尔伯特变换、相乘这三种方法去解调一个chirp信号（振幅调制了一下，用来判断幅度噪声的影响），这个chirp信号频率范围就是sstv的范围。
用python3运行，可以对比3种算法的效果。
import matplotlib.pyplot as plt import numpy as np from scipy.signal import hilbert, chirp import math from math import sin, cos, pi from scipy.signal.windows import hann frame_size = 50 over_lap = 10 duration = 1.0 fs = 8000.0 samples = int(fs*duration) t = np.arange(samples) / fs signal = chirp(t, 1500.0, t[-1], 2300.0) #1200~2300Hz signal *= (1.0 + 0.5 * np.sin(2.0*np.pi*3.0*t) ) def zero_cr(signal_data, frame_size, over_lap): wlen = len(signal_data) step = frame_size - over_lap print (wlen, step, math.ceil(wlen/step)) frame_num = int(math.ceil(wlen/step)) zcr = np.zeros((frame_num, 1)) for i in range(frame_num): cur_frame = signal_data[np.arange(i*step, min(i*step + frame_size, wlen))] cur_frame = cur_frame - np.mean(cur_frame) zcr[i] = sum(cur_frame[0:-2:1] * cur_frame[1:-1:1]<=0) * fs / frame_size /1.5 return zcr FIRLPCoeffs = [ +0.0001082461, +0.0034041195, +0.0063570207, +0.0078081648, +0.0060550614, -0.0002142384, -0.0104500335, -0.0211855480, -0.0264527776, -0.0201269304, +0.0004419626, +0.0312014771, +0.0606261038, +0.0727491887, +0.0537028370, -0.0004362161, -0.0779387981, -0.1511168919, -0.1829049634, -0.1390189257, -0.0017097774, +0.2201896764, +0.4894395006, +0.7485289338, +0.9357596142, +1.0040320616, +0.9357596142, +0.7485289338, +0.4894395006, +0.2201896764, -0.0017097774, -0.1390189257, -0.1829049634, -0.1511168919, -0.0779387981, -0.0004362161, +0.0537028370, +0.0727491887, +0.0606261038, +0.0312014771, +0.0004419626, -0.0201269304, -0.0264527776, -0.0211855480, -0.0104500335, -0.0002142384, +0.0060550614, +0.0078081648, +0.0063570207, +0.0034041195, +0.0001082461] xvFIRLP1 = [0 for _ in range(len(FIRLPCoeffs))] def FIRLowPass1(sampleIn): for i in range(0, len(xvFIRLP1)-1): xvFIRLP1[i] = xvFIRLP1[i+1] xvFIRLP1[len(xvFIRLP1)-1] = sampleIn / 5.013665674 sum = 0 for i in range(0, len(xvFIRLP1)): sum += (FIRLPCoeffs[i] * xvFIRLP1[i]) return sum xvFIRLP2 = [0 for _ in range(len(FIRLPCoeffs))] def FIRLowPass2(sampleIn): for i in range(0, len(xvFIRLP2)-1): xvFIRLP2[i] = xvFIRLP2[i+1] xvFIRLP2[len(xvFIRLP2)-1] = sampleIn / 5.013665674 sum = 0 for i in range(0, len(xvFIRLP2)): sum += (FIRLPCoeffs[i] * xvFIRLP2[i]) return sum xvMA1 = [0 for _ in range(9)] yvMA1prev = 0 def noiseReductionFilter1(sampleIn): global yvMA1prev for i in range(0, len(xvMA1)-1): xvMA1[i] = xvMA1[i+1] xvMA1[len(xvMA1)-1] = sampleIn yvMA1prev = yvMA1prev + xvMA1[len(xvMA1)-1] - xvMA1[0] return yvMA1prev xvMA2 = [0 for _ in range(9)] yvMA2prev = 0 def noiseReductionFilter2(sampleIn): global yvMA2prev for i in range(0, len(xvMA2)-1): xvMA2[i] = xvMA2[i+1] xvMA2[len(xvMA2)-1] = sampleIn yvMA2prev = yvMA2prev + xvMA2[len(xvMA2)-1] - xvMA2[0] return yvMA2prev oscPhase = 0 realPartPrev = 0 imaginaryPartPrev = 0 def fmDemodulate(sample): global oscPhase, realPartPrev, imaginaryPartPrev #global realPart, imaginaryPart oscPhase = oscPhase + (2 * pi * 2000) / fs realPart = cos(oscPhase) * sample #realPart_array[i] = cos(oscPhase) imaginaryPart = sin(oscPhase) * sample #imaginaryPart_array[i] = imaginaryPart if (oscPhase >= 2 * pi): oscPhase = oscPhase - 2 * pi realPart = FIRLowPass1(realPart) imaginaryPart = FIRLowPass2(imaginaryPart) realPart = noiseReductionFilter1(realPart) imaginaryPart = noiseReductionFilter2(imaginaryPart) demod_result = (imaginaryPart*realPartPrev-realPart*imaginaryPartPrev)/(realPart*realPart+imaginaryPart*imaginaryPart) realPartPrev = realPart imaginaryPartPrev = imaginaryPart demod_result = demod_result - 0.2335 - 0.694 luminance = int(round((demod_result/0.617)*800, 0)) luminance = 800-luminance #if (luminance > 255): # luminance = 255 #if (luminance < 0): # luminance = 0 return luminance realPart_array = np.zeros((len(signal), 1)) imaginaryPart_array = np.zeros((len(signal), 1)) frequency_multiply = np.zeros((len(signal), 1)) for i in range(0, len(signal)): frequency_multiply[i] = fmDemodulate(signal[i]) #-- Generate the Hilbert transformed vector of the samples analytic_signal = hilbert(signal) amplitude_envelope = np.abs(analytic_signal) instantaneous_phase = np.unwrap(np.angle(analytic_signal)) instantaneous_frequency = np.diff(instantaneous_phase) / (2.0*np.pi) * fs zcr = zero_cr(signal, frame_size, over_lap) #-- Plot the frequency over time graph plt.plot(signal) plt.ylabel('signal') plt.show() plt.plot(instantaneous_frequency) plt.ylabel('instantaneous_frequency') plt.show() plt.plot(zcr) plt.ylabel('zcr') plt.show() plt.plot(frequency_multiply) plt.ylabel('frequency_multiply') plt.show()
可以看到相乘法的效果和希尔伯特差不多，基本就是一条直线。
接下来我用相乘法和希尔伯特变换，分别解调了一个8000采样率，robot8bw制式的wav文件。
https://github.com/MossFrog/SSTV-Decoder
这个文件是从上面链接下载的，我用audacity改了example.wav的采样率，改到8000Hz。
from scipy.io.wavfile import read import matplotlib.pyplot as plt import pyaudio import numpy as np from math import sin, cos, pi from PIL import Image from scipy.signal import hilbert #-- Initial Variables SRate = 8000 FIRLPCoeffs = [ +0.0001082461, +0.0034041195, +0.0063570207, +0.0078081648, +0.0060550614, -0.0002142384, -0.0104500335, -0.0211855480, -0.0264527776, -0.0201269304, +0.0004419626, +0.0312014771, +0.0606261038, +0.0727491887, +0.0537028370, -0.0004362161, -0.0779387981, -0.1511168919, -0.1829049634, -0.1390189257, -0.0017097774, +0.2201896764, +0.4894395006, +0.7485289338, +0.9357596142, +1.0040320616, +0.9357596142, +0.7485289338, +0.4894395006, +0.2201896764, -0.0017097774, -0.1390189257, -0.1829049634, -0.1511168919, -0.0779387981, -0.0004362161, +0.0537028370, +0.0727491887, +0.0606261038, +0.0312014771, +0.0004419626, -0.0201269304, -0.0264527776, -0.0211855480, -0.0104500335, -0.0002142384, +0.0060550614, +0.0078081648, +0.0063570207, +0.0034041195, +0.0001082461] xvFIRLP1 = [0 for _ in range(len(FIRLPCoeffs))] def FIRLowPass1(sampleIn): for i in range(0, len(xvFIRLP1)-1): xvFIRLP1[i] = xvFIRLP1[i+1] xvFIRLP1[len(xvFIRLP1)-1] = sampleIn / 5.013665674 sum = 0 for i in range(0, len(xvFIRLP1)): sum += (FIRLPCoeffs[i] * xvFIRLP1[i]) return sum xvFIRLP2 = [0 for _ in range(len(FIRLPCoeffs))] def FIRLowPass2(sampleIn): for i in range(0, len(xvFIRLP2)-1): xvFIRLP2[i] = xvFIRLP2[i+1] xvFIRLP2[len(xvFIRLP2)-1] = sampleIn / 5.013665674 sum = 0 for i in range(0, len(xvFIRLP2)): sum += (FIRLPCoeffs[i] * xvFIRLP2[i]) return sum xvMA1 = [0 for _ in range(9)] yvMA1prev = 0 def noiseReductionFilter1(sampleIn): global yvMA1prev for i in range(0, len(xvMA1)-1): xvMA1[i] = xvMA1[i+1] xvMA1[len(xvMA1)-1] = sampleIn yvMA1prev = yvMA1prev + xvMA1[len(xvMA1)-1] - xvMA1[0] return yvMA1prev xvMA2 = [0 for _ in range(9)] yvMA2prev = 0 def noiseReductionFilter2(sampleIn): global yvMA2prev for i in range(0, len(xvMA2)-1): xvMA2[i] = xvMA2[i+1] xvMA2[len(xvMA2)-1] = sampleIn yvMA2prev = yvMA2prev + xvMA2[len(xvMA2)-1] - xvMA2[0] return yvMA2prev oscPhase = 0 realPartPrev = 0 imaginaryPartPrev = 0 def fmDemodulate(sample): global oscPhase, realPartPrev, imaginaryPartPrev oscPhase = oscPhase + (2 * pi * 2000) / SRate realPart = cos(oscPhase) * sample imaginaryPart = sin(oscPhase) * sample if (oscPhase >= 2 * pi): oscPhase = oscPhase - 2 * pi realPart = FIRLowPass1(realPart) imaginaryPart = FIRLowPass2(imaginaryPart) realPart = noiseReductionFilter1(realPart) imaginaryPart = noiseReductionFilter2(imaginaryPart) demod_result = (imaginaryPart*realPartPrev-realPart*imaginaryPartPrev)/(realPart*realPart+imaginaryPart*imaginaryPart) realPartPrev = realPart imaginaryPartPrev = imaginaryPart demod_result = demod_result - 0.2335 - 0.694 luminance = int(round((demod_result/0.617)*800, 0)) luminance = 800-luminance return luminance #-- Create the lines array (image) img = Image.new( 'L', (140,140), "white") pixelData = img.load() #-- Copy all the samples into an array input_data = read("example8000.wav") audio = input_data[1] frequency_multiply = np.zeros((len(audio), 1)) for i in range(0, len(audio)): frequency_multiply[i] = fmDemodulate(audio[i]) #-- Generate the Hilbert transformed vector of the samples analytic_signal = hilbert(audio) amplitude_envelope = np.abs(analytic_signal) instantaneous_phase = np.unwrap(np.angle(analytic_signal)) frequency_hilbert = np.diff(instantaneous_phase) / (2.0*np.pi) * SRate i_frequency = frequency_multiply #i_frequency = frequency_hilbert #-- Plot the frequency over time graph plt.plot(frequency_multiply) plt.ylabel('frequency_multiply') plt.show() plt.plot(frequency_hilbert) plt.ylabel('frequency_hilbert') plt.show() lineCount = 0 sampleBuffer = 0 #-- Search for the line sync pulses (Should be approx above 190 samples below 1300hz for i in range(0, len(i_frequency)-490): #-- If the frequency is below 1300Hz start counting if(i_frequency[i] < 1300): sampleBuffer = sampleBuffer + 1 #-- If over 207 samples below the given frequency are found increment the line count if(sampleBuffer > 29): lineCount = lineCount + 1 sampleBuffer = 0 lineBuffer = 0 lineIndex = 0 #-- Grab the whole line data for j in range(i, i + 490): lineBuffer = lineBuffer + 1 if(lineBuffer > 4): lineBuffer = 0 #-- Map each individual data point into a pixel value if(i_frequency[j] < 1500): pixelData[lineIndex,lineCount] = 0 elif(i_frequency[j] > 2300): pixelData[lineIndex,lineCount] = 255 #-- Normalize the value from frequency to Grayscale and map the pixel else: pixelData[lineIndex,lineCount] = np.int(((i_frequency[j] - 1500.0)/800.0)*255.0) #-- Skip a line lineIndex = lineIndex + 1 #-- Safety check to prevent extra line detection if(i_frequency[i] > 1300): sampleBuffer = 0 print("-- Line Count --") print(lineCount) print("-- Sample Rate --") print(SRate) print("-- Number of Samples --") print(len(i_frequency)) print("") print("== Operation Complete... ==") #-- Display the resulting image img.show()运行后发现，两种算法都能成功解调出图像。
hilbert frequency:
multiply frequency:
photo hilbert
photo multiply
所以接下来我打算用相乘法去解调robot8，并且用c++实现，这样方便搬入portapack。
c++ decode.cpp
#include <stdio.h> #include <stdint.h> #include <stdlib.h> #include <string.h> #include <math.h> #include <complex.h> #include <time.h> #include <iostream> #include <opencv2/core/core.hpp> #include <opencv2/highgui/highgui.hpp> #include <opencv2/imgproc/imgproc.hpp> #include <sys/stat.h> #include <fcntl.h> #include <sys/mman.h> using namespace cv; using namespace std; struct pcm *pcm; int SAMPLERATE = 8000; struct pcm { int (*rw)(struct pcm *, short *, int); void *data; }; struct wav_head { uint32_t ChunkID; uint32_t ChunkSize; uint32_t Format; uint32_t Subchunk1ID; uint32_t Subchunk1Size; uint16_t AudioFormat; uint16_t NumChannels; uint32_t SampleRate; uint32_t ByteRate; uint16_t BlockAlign; uint16_t BitsPerSample; uint32_t Subchunk2ID; uint32_t Subchunk2Size; }; struct wav { struct pcm base; void *p; short *b; size_t size; int index; int frames; int r; int c; }; int mmap_file_ro(void **p, char *name, size_t *size) { *size = 0; int fd = open(name, O_RDONLY); if (fd == -1) { perror("open"); return 0; } struct stat sb; if (fstat(fd, &sb) == -1) { perror("fstat"); //close(fd); return 0; } if (!S_ISREG(sb.st_mode)) { fprintf(stderr, "%s not a file\n", name); //close(fd); return 0; } *p = mmap(0, sb.st_size, PROT_READ, MAP_SHARED, fd, 0); if (*p == MAP_FAILED) { perror("mmap"); //close(fd); return 0; } /*if (close(fd) == -1) { perror ("close"); return 0; }*/ *size = sb.st_size; return 1; } int read_wav(struct pcm *pcm, short *buff, int frames) { struct wav *wav = (struct wav *)(pcm->data); if ((wav->index + frames) > wav->frames) return 0; memcpy(buff, wav->b + wav->index * wav->c, sizeof(short) * frames * wav->c); wav->index += frames; return 1; } int open_wav_read(struct pcm **p, char *name) { struct wav *wav = (struct wav *)malloc(sizeof(struct wav)); wav->base.rw = read_wav; wav->base.data = (void *)wav; if (!mmap_file_ro(&wav->p, name, &wav->size)) { fprintf(stderr, "couldnt open wav file %s!\n", name); free(wav); return 0; } struct wav_head *head = (struct wav_head *)wav->p; wav->b = (short *)(wav->p + sizeof(struct wav_head)); if (head->ChunkID != 0x46464952 || head->Format != 0x45564157 || head->Subchunk1ID != 0x20746d66 || head->Subchunk1Size != 16 || head->AudioFormat != 1 || head->Subchunk2ID != 0x61746164) { fprintf(stderr, "unsupported WAV file!\n"); free(wav); return 0; } if (head->BitsPerSample != 16) { fprintf(stderr, "only 16bit WAV supported!\n"); free(wav); return 0; } wav->index = 0; wav->frames = head->Subchunk2Size / (sizeof(short) * head->NumChannels); wav->r = head->SampleRate; wav->c = head->NumChannels; *p = &(wav->base); return 1; } double round(double r) { return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5); } double FIRLPCoeffs[51] = { +0.0001082461, +0.0034041195, +0.0063570207, +0.0078081648, +0.0060550614, -0.0002142384, -0.0104500335, -0.0211855480, -0.0264527776, -0.0201269304, +0.0004419626, +0.0312014771, +0.0606261038, +0.0727491887, +0.0537028370, -0.0004362161, -0.0779387981, -0.1511168919, -0.1829049634, -0.1390189257, -0.0017097774, +0.2201896764, +0.4894395006, +0.7485289338, +0.9357596142, +1.0040320616, +0.9357596142, +0.7485289338, +0.4894395006, +0.2201896764, -0.0017097774, -0.1390189257, -0.1829049634, -0.1511168919, -0.0779387981, -0.0004362161, +0.0537028370, +0.0727491887, +0.0606261038, +0.0312014771, +0.0004419626, -0.0201269304, -0.0264527776, -0.0211855480, -0.0104500335, -0.0002142384, +0.0060550614, +0.0078081648, +0.0063570207, +0.0034041195, +0.0001082461, }; double xvFIRLP1[51]; double FIRLowPass1(double sampleIn) { for (int i = 0; i < 50; i++) xvFIRLP1[i] = xvFIRLP1[i+1]; xvFIRLP1[50] = sampleIn / 5.013665674; double sum = 0; for (int i = 0; i <= 50; i++) sum += (FIRLPCoeffs[i] * xvFIRLP1[i]); return sum; } double xvFIRLP2[51]; double FIRLowPass2(double sampleIn) { for (int i = 0; i < 50; i++) xvFIRLP2[i] = xvFIRLP2[i+1]; xvFIRLP2[50] = sampleIn / 5.013665674; double sum = 0; for (int i = 0; i <= 50; i++) sum += (FIRLPCoeffs[i] * xvFIRLP2[i]); return sum; } // moving average double xvMA1[9]; double yvMA1prev = 0; double noiseReductionFilter1(double sampleIn) { for (int i = 0; i < 8; i++) xvMA1[i] = xvMA1[i+1]; xvMA1[8] = sampleIn; yvMA1prev = yvMA1prev+xvMA1[8]-xvMA1[0]; return yvMA1prev; } double xvMA2[9]; double yvMA2prev = 0; double noiseReductionFilter2(double sampleIn) { for (int i = 0; i < 8; i++) xvMA2[i] = xvMA2[i+1]; xvMA2[8] = sampleIn; yvMA2prev = yvMA2prev+xvMA2[8]-xvMA2[0]; return yvMA2prev; } double oscPhase = 0; double realPartPrev = 0; double imaginaryPartPrev = 0; int fmDemodulateLuminance(double sample) { oscPhase += (2 * M_PI * 2000) / SAMPLERATE; double realPart = cos(oscPhase) * sample; double imaginaryPart = sin(oscPhase) * sample; if (oscPhase >= 2 * M_PI) oscPhase -= 2 * M_PI; realPart = FIRLowPass1(realPart); imaginaryPart = FIRLowPass2(imaginaryPart); realPart = noiseReductionFilter1(realPart); imaginaryPart = noiseReductionFilter2(imaginaryPart); sample = (imaginaryPart*realPartPrev-realPart*imaginaryPartPrev)/(realPart*realPart+imaginaryPart*imaginaryPart); realPartPrev = realPart; imaginaryPartPrev = imaginaryPart; sample = sample - 0.2335 - 0.694; int luminance = (int)round((sample/0.617)*800); luminance = 800-luminance; return luminance; } int main(int argc, char **argv) { int pixelData[140][140]; char *pcm_name; int nframes = 61282; int lineCount = 0; int sampleBuffer = 0; pcm_name = "example8000_new.wav"; Mat frame; frame = Mat::zeros(67, 98, CV_8UC1); if (!open_wav_read(&pcm, pcm_name)) return 1; double sampleDouble[nframes]; double frequency_multiply[nframes]; for (int i; i < nframes - 1470; i++) { static short *pcm_buff; pcm_buff = (short *)malloc(sizeof(short) * 1); if (!read_wav(pcm, pcm_buff, 1)) { free(pcm_buff); return 0; } sampleDouble[i] = ((double)pcm_buff[0]) / 32768.0; frequency_multiply[i] = fmDemodulateLuminance(sampleDouble[i]); cout << "i: " << i << " freq: " << frequency_multiply[i] << endl; if (frequency_multiply[i] < 1300) { sampleBuffer = sampleBuffer + 1; } if (sampleBuffer > 29) { cout << "new line" << lineCount << endl; lineCount = lineCount + 1; sampleBuffer = 0; int lineIndex = 0; for (int j = i - 490 - 29; j < i; j=j+4) { if (frequency_multiply[j] < 1500) { pixelData[lineIndex][lineCount] = 0; frame.at<uchar>(lineCount, lineIndex) = 0; } else if (frequency_multiply[j] > 2300) { pixelData[lineIndex][lineCount] = 255; frame.at<uchar>(lineCount, lineIndex) = 255; } else { pixelData[lineIndex][lineCount] = int(((frequency_multiply[j] - 1500.0)/800.0)*255.0); frame.at<uchar>(lineCount, lineIndex) = int(((frequency_multiply[j] - 1500.0)/800.0)*255.0); } lineIndex = lineIndex + 1; imshow("frame", frame); if (waitKey(5) == 'q') { } } } if(frequency_multiply[i] > 1300) { sampleBuffer = 0; } } imwrite("frame.jpg", frame); return 0; }我这个版本的c++实现中，已经可以边解调边画图了，之前的robot8解调代码是先把所有wav数据全部fm解调，解调完毕后才去解码画图，实时性不行，我这里合并了解调和解码的过程，一边收一边绘图。
要用上面的代码，要下面2个命令，下面第一行是用来转换wav格式的。
sox example8000.wav -b 16 -e signed-integer example8000_new.wav g++ decode.cpp -o decode -lm `pkg-config --cflags --libs opencv`
