Research and Implementation of Key Technologies of Satellite Signal Spectrum Monitoring System - Dong Xue
Implementation of Rate Signal Method
The main idea here is to pay attention to the details of index, remember to subtract 1.
Instead of using czt to refine the spectrum, zoomfft is used to refine the spectrum. This method is easier to implement in engineering. All spectrum transformations can be done only by moving arrays, without exponential calculation. Reduce the amount of calculation. The zoomfft principle can be referred to. Here
Question: Why the refined results are not correct and deviate from the estimated results is still unclear. It is possible that approximation in some places will lead to this result.
% Test Rate Signal Method clc clear all close all sr = 1000; fs = 8e3; N = sr; srcData = randsrc(1, N, [0:3]); modData = pskmod(srcData, 4); % Pulse shaping filter ipoint = fs / sr; span = 21; pointFir = span * ipoint; upData = upsample(modData, ipoint); sendPulser = rcosdesign(0.5, span, ipoint, 'sqrt'); sendData = conv(upData, sendPulser); % Form the transmitting signal, here mainly subtract the time delay caused by convolution. See Digital Signal Processing for details. sendData = sendData(pointFir/2 + 1 : end - pointFir/2); % Computation rate signal for i = 1:length(sendData) %Get the Rate Signal if(i<length(sendData)) res(i) = (norm(sendData(i+1)-sendData(i)))*fs; %norm(x)seek x The 2 norm end end % Find the maximum power length of 2 less than length n=nextpow2(length(res)); %Find out 2^n >= length(res) Of n value n=n-1; len=2^n; % Change in data length and resolution df = fs / len; res=res(1:len); freResponse = fft(res); ampResponse = abs(freResponse); % Display spectrum, after removing 0 frequency, the largest frequency component is the coarse estimation of symbol rate. plot((0:length(ampResponse)-1)*df, ampResponse); % Removal of 0 frequency a = ampResponse(1); ampResponse(1) = 0; [~, index] = max(ampResponse); ampResponse(1) = a; % Here's a detail. Attention matlab From 1 to index Start acquiring, and you'll miss 0. % So our index When turning frequency, remember-1 coarseSymbolRate = (index-1) * df % Thinning spectrum % Set the refinement factor, which needs to be set to 2 n pow m = 64; n0 = len / m; % Move the right spectrum to the right of the zero frequency n02 = n0 / 2; right = index + n02; tmp = right + 1 - index; assert(tmp == n02+1); freResponse(1:tmp) = freResponse(index:right); % Left spectrum to negative frequency left = index - n02 + 1; tmp = n02 - 1; assert(tmp == (index - left)); freResponse((length(freResponse) - tmp+1):end) = freResponse(left:(index-1)); freResponse(right+1:end-tmp) = 0; % ifft timeSignal = ifft(freResponse); % Down sampling downTimeSignal = downsample(timeSignal, m); downTimeSignalFreqResponse = fft(downTimeSignal, len); leftFreq = (left-1) * df; rightFreq = (right-1) * df; downTimeSignalFreqResponse = circshift(downTimeSignalFreqResponse, floor(len/2)); f = linspace(leftFreq, rightFreq, len); ampDownTimeSignalFreqResponse = abs(downTimeSignalFreqResponse); figure() plot(f, ampDownTimeSignalFreqResponse); [~, index] = max(ampDownTimeSignalFreqResponse); df = (rightFreq - leftFreq) / len; accurateSr = (index - 1) * df + leftFreq % The maximum point in the spectrum is the result.