%%%% 本程序思想:
%2015.7.6 BPSK循环平稳
% 初步认定为频域加窗的SSCA算法,运用信号与其复共扼频谱相乘的方式来表达循环谱。被处理信号为运用BPSK调制产生的码源序
%%%% 列。关于频率f和循环频率alfa的表达式计算值得关注(只与采样频率有关)。
%数据设置部分
clc,clear all;
Rb = 2; %码元速率
Ts = 1/Rb; %码元间隔
L = 2^7; %单位时间内信号抽样点数
num = 16; %一个段落中码元个数
N = num*L; %一个段落内总的抽样点数 (4096)
fc = 12; %载波频率
dt = Ts/L; %抽样时间间隔
df = 1/(dt*N); %频率间隔
Bs = df*N/2; %带宽 BPSK的带宽是基带符号速率的2倍,如果知道基带符号速率,BPSK带宽就确定了
fs = Bs*2;
f = df/2-Bs:df:Bs; %抽样频率
t = dt/2:dt:num*Ts; %抽样时间 t= 1/(4*64) : 1/(2*64) : 8;
a=sign(randn(1,num)); %码源序列 generate the initial signal (length is 16)
%%%%%%%%%%%% BPSK调制信号产生
%通过基带成型滤波器后的信号 function
for rr=1:num
I((rr-1)*L+1:rr*L)=a(rr); % number of I is 1024 .extend the length of signal 'a' to 64 times length.
end
%最终调制信号 the signal modulated by BPSK carrier
s0 = I.*cos(2*pi*fc*t); % the number of 't'and 's' are the same.(1024)
figure(1)
plot(t,s0,'r-');
axis([0 0.5 -1 1]); %look 35/1024 percent of the plot.
%%%%%%%%%%%%加性高斯白噪声信道%%%%%%%%%%%%%%%%
snrdb = 8;
Ec = sum(abs(s0).^2)./length(s0);
snr_ratio = 10.^(snrdb/10);
noise_var = sqrt(Ec/(2*snr_ratio));
[s] = add_whitenoise2(s0,noise_var);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M = 2; % ??????????????????????????
Fs=fs;
N = (M*Fs)/df; % the new value of 'N'is 2 times the value of old 'N'. (2048)
%%%%N = pow2 (nextpow2(N)); % windowing record for FFT (2048) % NEXTPOW2(N) returns the first P such that 2^P >= abs(N)
% X = POW2(Y) for each element of Y is 2 raised to the power Y
%设置FFT点数为1024
X = fft(s0,N); % fft of the truncated(被删节的) (or zero padded) time series (2048 columns)
X = fftshift(X);% shift components of fft (2048 columns)
Xc = conj(X); % precompute the complex conjugate(复共轭) vector (2048 columns)
S = zeros (N,N); % size of the Spectral Correlation Density matrix
f = zeros (N,N); % size of the frequency matrix;
alfa = zeros (N,N); % size of the cycle frequency matrix
F = Fs/(2*N); % precompute(预先计算的) constants - F = Fs/(2*N);
G = Fs/N; % precompute constants - G = Fs/N;
m = -M/2+1:M/2; % set frequency smoothing window index % m(0)=0,m(1)=1.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Algorithm:
for k = 1:N % fix k (2048 columns) % the function of fix is choose the integral value towards '0'.
% computes vectors of f and alfa,
% store frequency and cycle frequency data for given k.
k1 = 1:N;
f(k,k1) = F*(k+k1-1) - Fs/2; % Computes f values and shift them to center in zero (f = (K+L)/2N) [1]
alfa(k,k1) = G*(k-k1 + N-1) - Fs; % Computes alfa values and shift them to center in zero (alfa = (K-L)/N) [1]
for k1 = 1:N %fix k1 = J
%calculate X(K+m) & conj (X(J+m)) for arguments of X(1:N) only
B = max(1-k, 1-k1); % Largest min of 1 <= (K+m)| (J+m) <= N
A = min (N-k, N-k1); % Smallest max of 1 <= (K+m)| (J+m) <= N
n = m((m<=A) & (m>=B)); % fix the index out of range problem by truncating the window. % m(0)=0,m(1)=1.
if isempty(n)
S(k,k1) = 0;
else
p = k+n; q = k1+n;
Y = X(p).*Xc(q);
S(k,k1) = sum(Y); % the power spetral Correlation destiny function
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% haha=length(sign([real(Y)]))
S = abs(S./max(max(S)));% normalize output matrix
alfa
f
figure(2);
mesh(alfa, f, S); grid;
ylabel('Frequency (Hz)');xlabel('Cycle frequency (Hz)');
title('BPSK循环谱三维图');
figure(3);
plot(alfa,S,'r'); % 'alfa'is cycle frequency.'f'is frequency.'S'is Spetral Correlation Density Function.
axis([-100 100 0 1]);
ylabel('S');xlabel('Cycle frequency (Hz)');
figure(4);
plot(f,S,'r');
axis([-100 100 0 1]);
ylabel('S');xlabel('Frequency (Hz)');
- 1
- 2
前往页