CN109450829B - Method and apparatus for estimating code rate of digital modulation signal - Google Patents

Method and apparatus for estimating code rate of digital modulation signal Download PDF

Info

Publication number
CN109450829B
CN109450829B CN201811350243.2A CN201811350243A CN109450829B CN 109450829 B CN109450829 B CN 109450829B CN 201811350243 A CN201811350243 A CN 201811350243A CN 109450829 B CN109450829 B CN 109450829B
Authority
CN
China
Prior art keywords
frequency
sequence
instantaneous
code rate
signal
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201811350243.2A
Other languages
Chinese (zh)
Other versions
CN109450829A (en
Inventor
王玉桃
王胜
张昊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing Changfeng Space Electronics Technology Co Ltd
Original Assignee
Nanjing Changfeng Space Electronics Technology Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing Changfeng Space Electronics Technology Co Ltd filed Critical Nanjing Changfeng Space Electronics Technology Co Ltd
Priority to CN201811350243.2A priority Critical patent/CN109450829B/en
Publication of CN109450829A publication Critical patent/CN109450829A/en
Application granted granted Critical
Publication of CN109450829B publication Critical patent/CN109450829B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04L—TRANSMISSION OF DIGITAL INFORMATION, e.g. TELEGRAPHIC COMMUNICATION
    • H04L25/00—Baseband systems
    • H04L25/02—Details ; arrangements for supplying electrical power along data transmission lines
    • H04L25/0262—Arrangements for detecting the data rate of an incoming signal

Landscapes

  • Engineering & Computer Science (AREA)
  • Power Engineering (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Signal Processing (AREA)
  • Digital Transmission Methods That Use Modulated Carrier Waves (AREA)

Abstract

The invention discloses a method and a device for estimating the code rate of a digital modulation signal, wherein the method comprises the following steps: for at a preset sampling rate fsExtracting instantaneous amplitude A (l) and instantaneous frequency f (l) of the sampled signal; after Haar wavelet transform, respectively performing discrete Fourier transformPerforming a Reed transform; extracting pulse peak frequency positions and sequencing the obtained positions; performing difference on the sequence, and respectively recording the numerical value with the most occurrence times in the sequence and the corresponding ticket number; changing the scale factor of the wavelet transformation according to a preset condition, and repeating the steps until the transformation scale a of the wavelet transformation is larger than the data length; respectively determining the value with the most votes in all the values in the instantaneous amplitude A (l) and the instantaneous frequency f (l) and the corresponding number of the votes; and comparing the number of votes obtained by the two values to finally determine the code rate. The method does not depend on prior knowledge, and has high recognition accuracy under the condition of lower signal-to-noise ratio; the proposed 'ticket number statistics' strategy can effectively solve the problem of wavelet transformation scale selection and has self-adaptability.

Description

Method and apparatus for estimating code rate of digital modulation signal
Technical Field
The invention relates to a method and a device for estimating the code rate of a digital modulation signal, belonging to the technical field of communication.
Background
The purpose of communication signal modulation and identification is to judge the modulation mode of a signal according to less prior information for a section of intercepted communication signals. With the rapid development of communication technology, the systems and modulation patterns of communication signals become more complex and diversified, and the signal environment becomes denser, which makes the identification of communication signals more difficult. The communication signals include two major types of modulation signals: when a decision theory is utilized, a threshold for distinguishing the analog modulation signal from the digital modulation signal is often difficult to set. Because the analog modulation signal has no code rate, and the estimation result is an arbitrary irregular value, the method for estimating the code rate can effectively distinguish analog modulation from digital modulation, and the code rate is an important characteristic for distinguishing two modulation modes.
The code rate has important significance for the extraction and analysis of the digital modulation signal characteristics and the identification and demodulation, in the electronic reconnaissance, the digital modulation signal is subjected to orthogonal frequency conversion and low-pass filtering to obtain a tie signal, and in order to analyze and identify the signal, the code rate of the digital modulation signal needs to be known. It is therefore necessary to correctly estimate the code rate of the digitally modulated signal.
There are various methods for estimating the code rate of a digitally modulated signal, the core idea of which is to estimate the code rate based on the spacing between abrupt changes, based on the presence of abrupt changes between symbols of the modulated signal. For digitally modulated baseband signals, delay multiplication can be used to obtain the code rate of the signal. After delay multiplication, the bilateral power spectral density of the output signal mainly comprises three terms, wherein the first term is a direct current component, the second term is a code rate and higher harmonics, and the third term is a continuous spectrum. The code rate of the modulated signal may be estimated based on the second term. But the delay length must be greater than half the symbol width, which is difficult to guarantee in practical applications.
The paper modulation identification study of radio signals (the authors' kingdom) discloses the use of amplitude modulated waveforms p (n) carrying digital information for code rate estimation. For ASK signals, p (n) is the instantaneous amplitude of the signal; for PSK signals, p (n) is the instantaneous phase of the signal; for FSK signals, p (n) is the instantaneous frequency of the signal. The process of the method can be summarized as follows:
1. solving an L-order differential sequence from p (n), and solving an absolute value to obtain a p' (n) sequence, wherein L is 15;
2. calculating the mean m of the p' (n) sequencepAt 3mpDeciding the p' (n) sequence for a threshold;
3. when p' (n) > 3mpThe decision device outputs 1; otherwise, the judger outputs 0; obtaining the sequence of step positions l (n)
4. Taking the center of the step zone as a symbol conversion point to obtain a symbol conversion position sequence lc(n) in the pair lc(n) taking the first order difference to obtain a difference sequence lcd(N) of length Nd；
5. Let Nm＝min[lcd(n)]Calculating the sequence ni(n)＝round(lcd(n)/Nm) Obtaining an estimate of the average symbol spacing
Figure BDA0001864666820000021
6. The code rate estimate of the acquired signal is:
Figure BDA0001864666820000022
wherein f isbIs the rate of the incoming data.
The method has the disadvantages that the threshold in the code rate estimation method is valued through experience, the theoretical basis is lacked, and the value of the threshold has no adaptability to modulation signals of different parameters;
in addition, since wavelet transform has a good detection performance for transient signals, code rate estimation based on wavelet transform is another idea. The document [ VSSR method based differentiation of digital and analog modulation signals-yupeng ] describes the method in detail, and the code rate estimation method can be summarized as follows:
1. carrying out wavelet transformation on the collected signals s (t) to obtain A (t);
2. obtaining envelope values after wavelet transformation of A (t) to obtain B (t);
3. fourier transform is carried out on the B (t) to obtain C (f);
4. and (3) calculating the distance d between two adjacent peaks in the | C (f) |, wherein the ratio of the distance d to the number of points of the | C (f) | is the ratio of the symbol rate to the sampling rate of the signal, so that the symbol rate of the signal can be calculated.
On the basis of the method of wavelet transform, the selection of the scale of the wavelet transform is not explained, if the code rate of the acquired signal is completely unknown, the selection of the appropriate scale of the wavelet transform is difficult, thereby causing the estimation error of the code rate;
in addition, the existing code rate estimation method is not comprehensive enough to analyze various digital modulation signals, generally performs code rate estimation on signals of a certain specific modulation mode, lacks a uniform code rate estimation method and may rely on prior knowledge, and if the modulation mode (amplitude, phase and frequency), carrier frequency and carrier phase of signals and the like need to be known in advance, the parameters are difficult to be accurately obtained under actual conditions;
disclosure of Invention
The invention aims to overcome the defects of the existing code rate estimation technology, provides a digital modulation signal code rate estimation method based on wavelet transformation without prior knowledge, and has high self-adaptive degree.
In order to solve the technical problems, the invention adopts the following technical scheme:
the method for estimating the code rate of the digital modulation signal comprises the following steps:
step 1: for at a preset sampling rate fsExtracting instantaneous amplitude A (L) and instantaneous frequency f (L) from the signal sampled by the total sampling point number L;
step 2: respectively performing wavelet transformation on the instantaneous amplitude A (l) and the instantaneous frequency f (l) according to a set wavelet transformation scale, wherein Haar wavelets are selected as mother functions of the wavelet transformation to obtain CA(l) And Cf(l)；
And step 3: to CA(l) Performing discrete Fourier transform to obtain YA(fk) (ii) a Performing discrete Fourier transform on the frequency f (l) to obtain Yf(fk)；
And 4, step 4: extraction of YA(fk) And Yf(fk) And sequencing the obtained positions to respectively obtain a sequence xA(N1) and sequence xf(N2), N1 and N2 are positive integers;
and 5: for sequence xA(N1) obtaining Dx by differenceA(N1) sequence for DxAThe values in (N1) were counted: select DxA(N1) the maximum number of repetitions, denoted XAiThe corresponding number of repetitions is NAi(ii) a For sequence xf(N2) obtaining Dx by differencef(N2) sequence for DxfThe values in (N2) were counted: select Dxf(N2) the maximum number of repetitions, denoted XfjCorresponding number of repetitions of Nfj(ii) a Wherein the subscripts i and j are integers and have an initial value of zero;
step 6: changing a wavelet transformation scale factor M according to a preset condition, enabling i to be equal to i plus 1, enabling j to be equal to j plus 1, repeating the steps 2-5 until the transformation scale a of the wavelet transformation is larger than the data length, and recording that i and j at the moment are equal to M;
and 7: sorting out NA0,NA1,NA2,…NAMAnd recording the index ma as ma, ma being an integer from 0 to M, XAmaIs recorded as RSA，NAmaIs recorded as ticket number PA(ii) a Select Nf0,Nf1,Nf2,…NfMAnd record its index as mf, mf being an integer from 0 to M, XfmfIs recorded as RSf，NfmfIs recorded as ticket number Pf；
And 8: if PAGreater than PfIf the code rate RS is RSAAnd the sampling rate fsThe ratio of the product of (d) to the total number of sample points L; otherwise, the code rate RS is RSfAnd the sampling rate fsThe product of (d) and the total number of sample points L.
Further, extracting the instantaneous amplitude a (l) and the instantaneous frequency f (l) comprises:
for continuous wave modulation, the expression of the modulated signal is:
Figure BDA0001864666820000056
in the formula (f)cRepresenting the frequency of the carrier wave, a (t) representing the instantaneous amplitude of the modulated signal; f (t) represents frequency;
Figure BDA0001864666820000057
is the phase;
performing Hilbert transform on the signal to obtain an orthogonal component of the signal, wherein the expression is as follows:
Figure BDA0001864666820000051
the analytical expression z (t) for the signal s (t) is obtained as s (t) + jv (t).
Still further, the method for extracting the instantaneous amplitude A (l) comprises the following steps:
Figure BDA0001864666820000052
after sampling the signal, the instantaneous amplitude a (L) is obtained in discrete form, L1, 2, … L.
Further, the method for extracting the instantaneous frequency f (l) comprises the following steps: extracting instantaneous phase, wherein the expression of the instantaneous phase is as follows:
Figure BDA0001864666820000053
after sampling the signal, obtaining the instantaneous phase of discrete form
Figure BDA0001864666820000054
Sequence, 1,2, … L;
to remove the phase folding, a corrected phase sequence c (L) is first calculated, L1, 2, … L:
Figure BDA0001864666820000055
the phase without folding is:
Figure BDA0001864666820000061
carrier frequency linear phase removal: estimating linear phase components using a phasor programming method
Figure BDA0001864666820000062
And make an error
Figure BDA0001864666820000063
Minimum, find C1And C2Two constants, expressed as:
Figure BDA0001864666820000064
wherein
Figure BDA0001864666820000065
In the case of a non-linear phase component,
Figure BDA0001864666820000066
a phase with no folding;
instantaneous frequency extraction, expressed as
Figure BDA0001864666820000067
Wherein f issIs the sampling rate.
Further, the wavelet transformation scale factor is m, and m is an integer not less than zero, then the expression of the transformation scale a of the wavelet transformation is: a 2m。
Still further, the initial wavelet transform scale factor m is set to zero.
Further, the wavelet transformation scale factor m is equal to m plus 1 when the wavelet transformation scale factor is changed.
Further, the pulse peak frequency position is extracted as follows:
|Y(fk) If the frequency corresponding to the code rate and the frequency multiplication position of the code rate have pulse spikes, the following expression is satisfied:
|Y(fk)|＞|Y(fk-1)|&|Y(fk)|＞|Y(fk+1)| (9)，
this indicates at frequency fkThere is a spike of the pulse.
In another aspect, the present invention provides a code rate estimation apparatus for digitally modulated signals, comprising:
a sampling unit: for determining the sampling rate fsSampling the signal by the total sampling point number L;
instantaneous amplitude and instantaneous frequency extraction unit: extracting instantaneous amplitude A (l) and instantaneous frequency f (l) of the adopted signal;
a wavelet transform unit: the wavelet transformation is respectively carried out on the instantaneous amplitude A (l) and the instantaneous frequency f (l) according to the set wavelet transformation scale, the Haar wavelet is selected as the wavelet transformation mother function to obtain CA(l) And Cf(l)；
Discrete Fourier transform unit: for pair CA(l) Performing discrete Fourier transform to obtain YA(fk) (ii) a Performing discrete Fourier transform on the frequency f (l)Leaf transformation to obtain Yf(fk) (ii) a Pulse spike frequency position extraction unit: for extracting YA(fk) And Yf(fk) And sequencing the obtained positions to respectively obtain a sequence xA(N1) and sequence xf(N1), N1 and N2 are positive integers;
a difference unit: for sequence xA(N1) obtaining Dx by differenceA(N1) sequence and sequence xf(N2) obtaining Dx by differencef(N2) sequence;
a ticket number counting unit: for DxAThe values in (N1) were counted: select DxA(N1) the maximum number of repetitions, denoted XAmThe corresponding number of repetitions is NAm(ii) a And for DxfThe values in (N2) were counted: select Dxf(N2) the maximum number of repetitions, denoted XfmCorresponding number of repetitions of Nfm；
Where the lower subscript m is a variable that is different for each cycle;
wavelet transformation scale factor setting unit: the method comprises the steps of setting a wavelet transformation scale, changing a wavelet transformation scale factor according to preset conditions until the transformation scale a of the wavelet transformation is larger than the data length, and recording M equal to the wavelet transformation scale factor M at the moment, wherein M is an integer not smaller than zero;
code rate calculation unit: for sorting out NA0,NA1,NA2,…NAMAnd recording the index ma as ma, ma being an integer from 0 to M, XAmaIs recorded as RSA，NAmaIs recorded as ticket number PA；
Sorting out Nf0,Nf1,Nf2,…NfMAnd record its index as mf, mf being an integer from 0 to M, XfmfIs recorded as RSf，NfmfIs recorded as ticket number Pf；
If PAGreater than PfIf the code rate RS is RSAAnd the sampling rate fsThe ratio of the product of (d) to the total number of sample points L; otherwise, the code rate RS isRSfAnd the sampling rate fsThe product of (d) and the total number of sample points L.
The invention achieves the following beneficial effects:
1. the types of digital modulation signals are various, and the common signals are comprehensively analyzed and code rate estimated;
2. the code rate estimation is carried out under the condition that the signal is completely unknown without relying on prior knowledge, including a signal modulation mode, a signal-to-noise ratio, a carrier frequency, a carrier phase, code element synchronization and the like, and the identification accuracy is high under the condition of a lower signal-to-noise ratio;
3. the 'ticket number statistics' strategy provided by the invention can effectively solve the problem of wavelet transformation scale selection;
4. the method has self-adaptability, and can automatically select a proper scale for different code rates of a digital modulation signal; the algorithm rule is simple and convergent, the estimation of the code rate can be completed in a short time, and the method is suitable for online analysis.
Drawings
FIG. 1 is a schematic flow diagram of the process of the present invention.
Detailed Description
The invention is further described below with reference to the accompanying drawings. The following examples are only for illustrating the technical solutions of the present invention more clearly, and the protection scope of the present invention is not limited thereby.
Fig. 1 is a flow chart illustrating the method of the present invention, fig. 1 shows the method of code rate estimation of a digital modulation signal of the present invention,
step 1: for at a preset sampling rate fsExtracting instantaneous amplitude A (L) and instantaneous frequency f (L) from the signal after sampling by the total sampling point number L, wherein L is a discrete independent variable after sampling; because the symbol cannot be determined to be modulated in amplitude or frequency without knowing the signal, the method of the invention processes in amplitude and frequency:
the specific method comprises the following steps:
for continuous wave modulation, the general expression of a modulated signal is:
Figure BDA0001864666820000091
in the formula (f)cRepresenting the frequency of the carrier. The modulated signal may be "parasitic" in the instantaneous amplitude a (t), frequency f (t), and phase of the modulated signal, respectively
Figure BDA0001864666820000092
In (1).
(1) And performing Hilbert transform on the signal to obtain orthogonal components of the signal.
Figure BDA0001864666820000093
The hubert transform actually shifts the original signal by pi/2 phase, making it an own orthogonal pair. The analytical expression z (t) of the signal s (t) is s (t) + jv (t).
(2) Instantaneous amplitude extraction
Signal envelope instantaneous amplitude estimation:
Figure BDA0001864666820000094
(3) instantaneous phase extraction
Instantaneous phase estimation:
Figure BDA0001864666820000101
(4) after sampling the signal, the instantaneous amplitude A (l) sequence and instantaneous phase in discrete form are obtained
Figure BDA0001864666820000102
Sequence, L ═ 1,2, … L
(5) Dephasing folding
Since the actual phase is calculated modulo 2 pi, after the instantaneous phase is determined, the phase needs to be unfolded and the unfolded phase restored, for which purpose a correction phase sequence c (L) is first calculated, L being 1,2, … L:
Figure BDA0001864666820000103
the phase without folding is:
Figure BDA0001864666820000104
(6) linear phase of frequency-carrier removal
The phase of the signal comprising two parts, the carrier-frequency-induced linear phase component
Figure BDA0001864666820000105
And non-linear phase component of modulation phase
Figure BDA0001864666820000106
Thus the non-linear phase component
Figure BDA0001864666820000107
Can be calculated from the following formula:
Figure BDA0001864666820000108
at carrier frequency fcIn case of not being exactly known, the linear phase component can be estimated by using a phase programming method
Figure BDA0001864666820000109
And make an error
Figure BDA00018646668200001010
At minimum, C can be found1And C2Two constants, so there are:
Figure BDA00018646668200001011
(7) instantaneous frequency extraction
The instantaneous frequency can be calculated from the derivative of the instantaneous phase, and the digital signal can be directly differentiated, that is:
Figure BDA0001864666820000111
step 2: respectively performing wavelet transformation on the instantaneous amplitude A (l) and the instantaneous frequency f (l) according to a set wavelet transformation scale, wherein the Haar wavelet is simple, fast to calculate and most common, so the Haar wavelet is selected as the wavelet transformation mother function to obtain CA(l) And Cf(l)；
The scale factor of the wavelet transformation is m, and m is an integer not less than zero, then the expression of the transformation scale a of the wavelet transformation is: a 2m。
It can be seen from fig. 1 that the wavelet transform scale factor is set to zero when setting the initial wavelet transform scale.
And step 3: to CA(l) Performing discrete Fourier transform to obtain YA(fk) (ii) a Performing discrete Fourier transform on the frequency f (l) to obtain Yf(fk)；
And 4, step 4: extraction of YA(fk) And Yf(fk) And sequencing the obtained positions to respectively obtain a sequence xA(N1) and sequence xf(N2), N1 and N2 are positive integers;
YA(fk) And Yf(fk) A pulse peak appears at the frequency corresponding to the code rate and the code rate frequency multiplication position, and the frequency position corresponding to the peak is recorded; the resulting positions may be sorted in both ascending and descending order.
And 5: for sequence xA(N1) obtaining Dx by differenceA(N1) sequence for DxAThe values in (N1) were counted: select DxA(N1) the maximum number of repetitions, denoted XAiThe corresponding number of repetitions is NAi(ii) a For sequence xf(N2) obtaining Dx by differencef(N2) sequence for DxfThe values in (N2) were counted: select Dxf(N2) medium weightThe value with the most repetition number is recorded as XfjCorresponding number of repetitions of Nfj(ii) a Wherein the subscripts i and j are integers and have an initial value of zero;
step 6: changing a wavelet transformation scale factor M according to a preset condition, enabling i to be equal to i plus 1, enabling j to be equal to j plus 1, repeating the steps 2-5 until the transformation scale a of the wavelet transformation is larger than the data length, and recording that i and j at the moment are equal to M;
preferably, the wavelet transform scale factor is changed such that the wavelet transform scale factor m is equal to m plus 1.
And 7: sorting out NA0,NA1,NA2,…NAMAnd recording the index ma as ma, ma being an integer from 0 to M, XAmaIs recorded as RSA，NAmaIs recorded as ticket number PA(ii) a Sorting out Nf0,Nf1,Nf2,…NfMAnd record its index as mf, mf being an integer from 0 to M, XfmfIs recorded as RSf，NfmfIs recorded as ticket number Pf；
And 8: if PAGreater than PfIf the code rate RS is RSAAnd the sampling rate fsThe ratio of the product of (d) to the total number of sample points L; otherwise, the code rate RS is RSfAnd the sampling rate fsThe product of (d) and the total number of sample points L.
Principle of wavelet transform
Defining: psi (t) epsilon L2(R) Fourier transform thereof
Figure BDA0001864666820000121
When in use
Figure BDA0001864666820000122
The allowance condition (complete reconstruction condition or identity resolution condition) is satisfied:
Figure BDA0001864666820000124
let ψ (t) be a basic wavelet or mother wavelet. Translating and stretching the mother function to obtain:
Figure BDA0001864666820000123
it is referred to as a wavelet series. Where a is the scaling factor and b is the translation factor. For an arbitrary function x (t) e L2(R) continuous wavelet transform:
Figure BDA0001864666820000131
wavelet transform of the corresponding sampled signal:
Figure BDA0001864666820000132
in the present invention, x (l) corresponds to the instantaneous amplitude A (l) and the instantaneous frequency f (l).
Since the base wavelet acts as an observation window for the analyzed signal in the wavelet transform, ψ (t) should also satisfy the constraint condition of the general function:
Figure BDA0001864666820000133
generally, the scale parameter a and the translation parameter b in the continuous wavelet transform are respectively taken as discrete forms:
Figure BDA0001864666820000134
j ∈ Z is called discrete wavelet transform.
(2) Haar mother wavelet
The mathematical expression of the Haar mother wavelet is
Figure BDA0001864666820000135
Extraction method of pulse peak frequency position
|Y(fk) If the pulse peak appears at the frequency corresponding to the code rate and the frequency multiplication position of the code rate, | the following rule is satisfied:
|Y(fk)|＞|Y(fk-1)|&|Y(fk)|＞|Y(fk+1)| (17)
this indicates at frequency fkIs in existence of pulse spike'&'means AND gate'. To eliminate the effects of noise fluctuations, only the 3dB bandwidth is considered for the spike in the frequency in the algorithm implementation.
In another aspect, the present invention provides a digitally modulated signal code rate estimation apparatus, comprising:
a sampling unit: for determining the sampling rate fsSampling the signal by the total sampling point number L;
instantaneous amplitude and instantaneous frequency extraction unit: extracting instantaneous amplitude A (l) and instantaneous frequency f (l) of the adopted signal;
a wavelet transform unit: the wavelet transformation is respectively carried out on the instantaneous amplitude A (l) and the instantaneous frequency f (l) according to the set wavelet transformation scale, the Haar wavelet is selected as the wavelet transformation mother function to obtain CA(l) And Cf(l)；
Discrete Fourier transform unit: for pair CA(l) Performing discrete Fourier transform to obtain YA(fk) (ii) a Performing discrete Fourier transform on the frequency f (l) to obtain Yf(fk) (ii) a Pulse spike frequency position extraction unit: for extracting YA(fk) And Yf(fk) And sequencing the obtained positions to respectively obtain a sequence xA(N1) and sequence xf(N1), N1 and N2 are positive integers;
a difference unit: for sequence xA(N1) obtaining Dx by differenceA(N1) sequence and sequence xf(N2) obtaining Dx by differencef(N2) sequence;
a ticket number counting unit: for DxAThe values in (N1) were counted: select DxA(N1) the maximum number of repetitions, denoted XAmThe corresponding number of repetitions is NAm(ii) a And for DxfThe values in (N2) were counted: select Dxf(N2) the maximum number of repetitions, denoted XfmCorresponding number of repetitions of Nfm；
Where the lower subscript m is a variable that is different for each cycle;
wavelet transformation scale factor setting unit: the method comprises the steps of setting a wavelet transformation scale, changing a wavelet transformation scale factor according to preset conditions until the transformation scale a of the wavelet transformation is larger than the data length, and recording M equal to the wavelet transformation scale factor M at the moment, wherein M is an integer not smaller than zero;
code rate calculation unit: for sorting out NA0,NA1,NA2,…NAMAnd recording the index ma as ma, ma being an integer from 0 to M, XAmaIs recorded as RSA，NAmaIs recorded as ticket number PA；
Sorting out Nf0,Nf1,Nf2,…NfMAnd record its index as mf, mf being an integer from 0 to M, XfmfIs recorded as RSf，NfmfIs recorded as ticket number Pf；
If PAGreater than PfIf the code rate RS is RSAAnd the sampling rate fsThe ratio of the product of (d) to the total number of sample points L; otherwise, the code rate RS is RSfAnd the sampling rate fsThe product of (d) and the total number of sample points L.
Further, the wavelet transform unit sets the initial wavelet transform scale factor m to zero, and the expression of the transform scale a of the wavelet transform is: a 2m。
In the specific embodiment, the digital modulation signals 2ASK,4ASK, BPSK, QPSK,8PSK, OQPSK, PI/4-DQPSK,2FSK,4FSK,16QAM, and 64QAM, which are commonly used in the actual communication system, are taken as examples, and the code rate is automatically estimated by computer simulation, and a random sequence is used as the source of the digital modulation signals. The sequence of additive white gaussian noise AWGN is generated using a normally distributed random number generator with a mean of 0. The signal-to-noise ratio SNR varies from 5dB to 15dB, stepping by 1 dB. 20000 points are used for each identification sample, and each identification sample is simulated 400 times. The simulation platform operating system is WIN7, programmed with software MATLAB. Simulation shows that when the signal-to-noise ratio reaches 6dB, the estimation accuracy of the code rate is up to 99%.
Digital modulation technology is most commonly used in modern communication systems, and therefore, identification of a digital modulation mode and estimation of a key parameter code rate play a crucial role in analysis and feature extraction of communication signals. The wavelet transform-based digital modulation signal code rate estimation method does not need any prior knowledge, can self-adaptively determine the scale of wavelet transform according to a voting strategy, and has simple and convergent algorithm rules. The simulation result of common digital modulation signals shows that the estimation accuracy of the code rate can still reach more than 99 percent under the condition of low signal-to-noise ratio
The above description is only a preferred embodiment of the present invention, and it should be noted that, for those skilled in the art, several modifications and variations can be made without departing from the technical principle of the present invention, and these modifications and variations should also be regarded as the protection scope of the present invention.

Claims (10)

1. The method for estimating the code rate of the digital modulation signal is characterized by comprising the following steps:
step 1: for at a preset sampling rate fsExtracting instantaneous amplitude A (L) and instantaneous frequency f (L) from the signal sampled by the total sampling point number L;
step 2: respectively performing wavelet transformation on the instantaneous amplitude A (l) and the instantaneous frequency f (l) according to a set wavelet transformation scale, wherein Haar wavelets are selected as mother functions of the wavelet transformation to obtain CA(l) And Cf(l)；
And step 3: to CA(l) Performing discrete Fourier transform to obtain YA(fk) (ii) a Performing discrete Fourier transform on the frequency f (l) to obtain Yf(fk)；
And 4, step 4: extraction of YA(fk) And Yf(fk) And sequencing the obtained positions to respectively obtain a sequence xA(N1) and sequence xf(N2), N1 and N2 are positive integers;
and 5: for sequence xA(N1) obtaining Dx by differenceA(N1) sequence for DxAThe values in (N1) were counted: select DxA(N1) the maximum number of repetitions, denoted XAiThe corresponding number of repetitions is NAi(ii) a For sequence xf(N2) obtaining Dx by differencef(N2) sequence for DxfThe values in (N2) were counted: select Dxf(N2) the maximum number of repetitions, denoted XfjCorresponding number of repetitions of Nfj(ii) a Wherein the subscripts i and j are integers and have an initial value of zero;
step 6: changing a wavelet transformation scale factor M according to a preset condition, enabling i to be equal to i plus 1, enabling j to be equal to j plus 1, repeating the steps 2-5 until the transformation scale a of the wavelet transformation is larger than the data length, and recording that i and j at the moment are equal to M;
and 7: sorting out NA0,NA1,NA2,…NAMAnd recording the index ma as ma, ma being an integer from 0 to M, XAmaIs recorded as RSA，NAmaIs recorded as ticket number PA(ii) a Sorting out Nf0,Nf1,Nf2,…NfMAnd record its index as mf, mf being an integer from 0 to M, XfmfIs recorded as RSf，NfmfIs recorded as ticket number Pf；
And 8: if PAGreater than PfIf the code rate RS is RSAAnd the sampling rate fsThe ratio of the product of (d) to the total number of sample points L; otherwise, the code rate RS is RSfAnd the sampling rate fsThe product of (d) and the total number of sample points L.
2. The method of claim 1 wherein extracting the instantaneous amplitude a (l) and the instantaneous frequency f (l) comprises:
for continuous wave modulation, the expression of the modulated signal is:
Figure FDA0001864666810000021
in the formula (f)cRepresenting the frequency of the carrier wave, a (t) representing the instantaneous amplitude of the modulated signal; f (t) represents frequency;
Figure FDA0001864666810000022
is the phase;
performing Hilbert transform on the signal to obtain an orthogonal component of the signal, wherein the expression is as follows:
Figure FDA0001864666810000023
the analytical expression z (t) for the signal s (t) is obtained as s (t) + jv (t).
3. The method of claim 2 wherein the instantaneous amplitude a (l) is extracted by:
Figure FDA0001864666810000024
after sampling the signal, the instantaneous amplitude a (L) is obtained in discrete form, L1, 2, … L.
4. The method of claim 2 wherein the method of extracting the instantaneous frequency f (l) comprises: extracting instantaneous phase, wherein the expression of the instantaneous phase is as follows:
Figure FDA0001864666810000025
after sampling the signal, obtaining the instantaneous phase of discrete form
Figure FDA0001864666810000026
Sequence of，l＝1,2,…L；
To remove the phase folding, a corrected phase sequence c (L) is first calculated, L1, 2, … L:
Figure FDA0001864666810000031
the phase without folding is:
Figure FDA0001864666810000032
carrier frequency linear phase removal: estimating linear phase components using a phasor programming method
Figure FDA0001864666810000033
And make an error
Figure FDA0001864666810000034
Minimum, find C1And C2Two constants, expressed as:
Figure FDA0001864666810000035
wherein
Figure FDA0001864666810000036
In the case of a non-linear phase component,
Figure FDA0001864666810000037
a phase with no folding;
instantaneous frequency extraction, expressed as
Figure FDA0001864666810000038
Wherein f issIs the sampling rate.
5. The method according to claim 1, wherein the wavelet transform scale factor is m, and m is an integer not less than zero, and the transform scale a of the wavelet transform is expressed as: a 2m。
6. The method according to claim 5, wherein the initial wavelet transform scale factor m is set to zero.
7. The method of claim 5, wherein the wavelet transform scale factor m is equal to m plus 1 when the wavelet transform scale factor is changed.
8. The method of claim 1, wherein the pulse peak frequency location is extracted by:
|Y(fk) If the frequency corresponding to the code rate and the frequency multiplication position of the code rate have pulse spikes, the following expression is satisfied:
|Y(fk)|＞|Y(fk-1)|&|Y(fk)|＞|Y(fk+1)| (9)，
this indicates at frequency fkThere is a spike of the pulse.
9. A digitally modulated signal code rate estimation apparatus, comprising:
a sampling unit: for determining the sampling rate fsSampling the signal by the total sampling point number L;
instantaneous amplitude and instantaneous frequency extraction unit: extracting instantaneous amplitude A (l) and instantaneous frequency f (l) of the adopted signal;
a wavelet transform unit: the wavelet transformation is respectively carried out on the instantaneous amplitude A (l) and the instantaneous frequency f (l) according to the set wavelet transformation scale, the Haar wavelet is selected as the wavelet transformation mother function to obtain CA(l) And Cf(l)；
Discrete Fourier transform unit: for pair CA(l) Performing discrete Fourier transform to obtain YA(fk) (ii) a Performing discrete Fourier transform on the frequency f (l) to obtain Yf(fk) (ii) a Pulse spike frequency position extraction unit: for extracting YA(fk) And Yf(fk) And sequencing the obtained positions to respectively obtain a sequence xA(N1) and sequence xf(N1), N1 and N2 are positive integers;
a difference unit: for sequence xA(N1) obtaining Dx by differenceA(N1) sequence and sequence xf(N2) obtaining Dx by differencef(N2) sequence;
a ticket number counting unit: for DxAThe values in (N1) were counted: select DxA(N1) the maximum number of repetitions, denoted XAmThe corresponding number of repetitions is NAm(ii) a And for DxfThe values in (N2) were counted: select Dxf(N2) the maximum number of repetitions, denoted XfmCorresponding number of repetitions of Nfm；
Where the lower subscript m is a variable that is different for each cycle;
wavelet transformation scale factor setting unit: the method comprises the steps of setting a wavelet transformation scale, changing a wavelet transformation scale factor according to preset conditions until the transformation scale a of the wavelet transformation is larger than the data length, and recording M equal to the wavelet transformation scale factor M at the moment, wherein M is an integer not smaller than zero;
code rate calculation unit: for sorting out NA0,NA1,NA2,…NAMAnd recording the index ma as ma, ma being an integer from 0 to M, XAmaIs recorded as RSA，NAmaIs recorded as ticket number PA；
Sorting out Nf0,Nf1,Nf2,…NfMAnd record its index as mf, mf being an integer from 0 to M, XfmfIs recorded as RSf，NfmfIs recorded as ticket number Pf；
If PAGreater than PfIf the code rate RS is RSAAnd the sampling rate fsThe ratio of the product of (d) to the total number of sample points L; otherwise, the code rate RS is RSfAnd the sampling rate fsThe product of (d) and the total number of sample points L.
10. The apparatus for code rate estimation of digitally modulated signals according to claim 9, comprising: the wavelet transformation unit sets an initial wavelet transformation scale factor m to be zero, and the expression of the transformation scale a of the wavelet transformation is as follows: a 2m。
CN201811350243.2A 2018-11-14 2018-11-14 Method and apparatus for estimating code rate of digital modulation signal Active CN109450829B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811350243.2A CN109450829B (en) 2018-11-14 2018-11-14 Method and apparatus for estimating code rate of digital modulation signal

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811350243.2A CN109450829B (en) 2018-11-14 2018-11-14 Method and apparatus for estimating code rate of digital modulation signal

Publications (2)

Publication Number Publication Date
CN109450829A CN109450829A (en) 2019-03-08
CN109450829B true CN109450829B (en) 2020-12-25

Family

ID=65551756

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811350243.2A Active CN109450829B (en) 2018-11-14 2018-11-14 Method and apparatus for estimating code rate of digital modulation signal

Country Status (1)

Country Link
CN (1) CN109450829B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110690931B (en) * 2019-10-11 2020-12-22 中国科学院软件研究所 Digital signal adaptive code rate estimation method and device based on multi-wavelet-base combination
CN113259027B (en) * 2021-05-12 2022-09-06 深圳华创电科技术有限公司 Method for calculating hostile-my identification signal code rate based on Haar transform
CN116032709B (en) * 2022-12-06 2024-04-12 中国电子科技集团公司第三十研究所 Method and device for blind demodulation and modulation feature analysis of FSK signal without priori knowledge

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103441964B (en) * 2013-08-30 2016-05-18 清华大学 A kind of communication signal parameters method of estimation and device based on accidental resonance
CN103490824A (en) * 2013-09-04 2014-01-01 中国电子科技集团公司第四十一研究所 Reference signal recovery method for EVM analysis
CN103841600B (en) * 2014-03-13 2017-02-22 中国人民解放军信息工程大学 Signal element rate estimation method and electromagnetic spectrum monitoring system of sensor network
KR101611534B1 (en) * 2015-08-28 2016-04-11 한화탈레스 주식회사 Method for symbol rate estimation
CN106209701B (en) * 2016-07-01 2019-08-02 中国人民解放军信息工程大学 MFSK signal code rate-estimation method and device under Alpha Stable distritation noise circumstance
CN108270700A (en) * 2016-12-30 2018-07-10 中国航天科工集团八五研究所 A kind of improved digital signal symbol rate feature extraction algorithm
CN108282428A (en) * 2018-01-23 2018-07-13 大连理工大学 A kind of automatic identifying method of common modulation mode of communication signal without priori

Also Published As

Publication number Publication date
CN109450829A (en) 2019-03-08

Similar Documents

Publication Publication Date Title
Wong et al. Specific emitter identification using convolutional neural network-based IQ imbalance estimators
CN109450829B (en) Method and apparatus for estimating code rate of digital modulation signal
CN107948107B (en) Digital modulation signal classification method based on joint features
CN110113278B (en) Modulation mode identification method based on all-digital receiver
CN110690931B (en) Digital signal adaptive code rate estimation method and device based on multi-wavelet-base combination
CN111935046B (en) Low-complexity frequency shift keying signal symbol rate estimation method
CN113037663B (en) Improved code element rate estimation method suitable for non-constant envelope signal
CN100521670C (en) Detecting and analyzing method for multi system frequency shift key control signal
CN104363194A (en) PSK (phase shift keying) modulation recognition method based on wave form transformation
CN105785324A (en) MGCSTFT-based chirp signal parameter estimation method
CN116257752A (en) Signal modulation pattern recognition method
CN105429719A (en) Strong interference signal detection method based on power spectrum and multiple dimensioned wavelet transformation analysis
Majhi et al. Novel blind modulation classification of circular and linearly modulated signals using cyclic cumulants
Wong et al. Emitter identification using CNN IQ imbalance estimators
CN108737302A (en) The symbol rate estimation method and its device of accidental resonance joint wavelet transformation under Low SNR
Haq et al. Recognition of digital modulated signals based on statistical parameters
Ghauri KNN based classification of digital modulated signals
Hassanpour et al. A robust algorithm based on wavelet transform for recognition of binary digital modulations
CN116032709B (en) Method and device for blind demodulation and modulation feature analysis of FSK signal without priori knowledge
CN109660475B (en) A kind of non-cooperation phase code water sound communication signal autonomous identifying method
CN116506273A (en) Novel MPSK modulation signal identification and classification method
CN109525528B (en) Image domain signal identification method facing MQAM modulation signal
CN115664905A (en) Wi-Fi equipment identification system and method based on multi-domain physical layer fingerprint characteristics
Wang et al. Blind symbol rate estimation of satellite communication signal by Haar wavelet transform
CN115460048A (en) MSK modulation identification method, medium and device based on code element rate

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant