扩频调制matlab仿真

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
%%
%1.测试产生m序列，并且查看其特性
clear;
clc;
close all;
%通过异或的方式产生m序列，这里选择了15位单极性m序列
M=50;
m_sequ=[];
for i=1:M
    register = m_sequence([1 0 0 1]);
    %产生m序列数组，用m_sequ保存起来
    m_sequ = [m_sequ register];
end
%将单极性的m序列转换成为双极性的m序列
bipolar_m_sequ = 1-2*m_sequ;
num=1:length(bipolar_m_sequ);

figure(1)
subplot(2,1,1)
stem(num-1,bipolar_m_sequ);
axis([0,15,-1.2,1.2]);%展示15位M序列
xlabel('num');
ylabel('M序列');
title('双极性15位M序列');
sample_num=4096;%设定进行fft变换的点数为4096个

subplot(2,1,2)
m_sequ_fft=fft(bipolar_m_sequ, sample_num);%双极性M序列进行fft变换
magnitude=abs(m_sequ_fft);%采样后的振幅magnitude
frequency_sequence=(1:sample_num)*200/sample_num;%频率序列，按采样频率为200进行采样
%进行fft变换得到的幅值不是真实的，需要将变换后的结果乘2再除以个数N
plot(frequency_sequence,magnitude(1:sample_num)*2/sample_num);
xlabel('Hz');
axis([100,115,0,0.2]);
title('绘出M序列随频率变化的振幅');


%%
%2.生成随机二进制比特序列，并进行扩频编码
N=50;
%生成50位随机二进制比特流
random_bitlist=ceil(rand(1,N)-0.5);
y=zeros(1,15*N);%生成频率信息码
for i=1:N
   k=15*i-14;
   y(k:k+15)=random_bitlist(i);
end
m = 750;%进行m序列编码（异或操作）
for i=1:m
    m_modulation_bitlist(i)=double(xor(m_sequ(i),y(i)));
end

figure(2)
subplot(2,1,1)
t=0:N-1;
stem(t,random_bitlist);
axis([0,50,-0.5,1.5])
title('扩频前二进制比特流');

subplot(2,1,2)
t=1:m;
stem(t,m_modulation_bitlist);
axis([0,200,-0.5,1.5]);
title('扩频后的二进制比特流');

%%
%3.对扩频前后信号进行BPSK调制，观察其时域波形

%设定采样频率
ts=0:1e-5:7.5-0.00001;
%每个符号复制1000次，将本来是冲激信号的m序列调制后比特流转变成矩形信号
m_modulation_bitlist_1=rectpulse(m_modulation_bitlist,1000);
%扩频后信号BPSK调制时域波形
fs=2000;%设定bpsk调制频率为2000hz
m_modulation_bitlist_bpsk=(1-2.*m_modulation_bitlist_1).*cos(2*pi*fs*ts);

m_modulation_bitlist_2=rectpulse(random_bitlist,15000);
%无扩频信号BPSK调制时域波形
m_modulation_bitlis_unbpsk=(1-2.*m_modulation_bitlist_2).*cos(2*pi*fs*ts);

figure(3)
subplot(2,1,2)
plot(ts,m_modulation_bitlist_bpsk);
xlabel('t/s');
axis([0, 0.05, -1.2, 1.2]);
title('扩频后bpsk信号时域波形');

subplot(2,1,1)
plot(ts,m_modulation_bitlis_unbpsk);
xlabel('s');
axis([0, 0.05, -1.2, 1.2]);
title('扩频前bpsk信号时域波形')



%%
%4.计算并观察扩频前后BPSK调制信号的频谱
%对信号采用800000点fft计算，得到频谱
figure(4)
N=800000;

y_unbpsk=fft(m_modulation_bitlis_unbpsk,N);
magnitude_unbpsk=abs(y_unbpsk);
frequency_unbpsk=(1:N/2)*100000/N;

y_bpsk=fft(m_modulation_bitlist_bpsk, N);
magnitude_bpsk=abs(y_bpsk);
frequency_bpsk=(1:N/2)*100000/N;

subplot(2,1,1)
plot(frequency_unbpsk,magnitude_unbpsk(1:N/2)*2/N);
%扩频前调制信号频率在2000hz
xlabel('Hz');
axis([1700,2300,0,0.8]);
title('扩频前调制信号频谱');

subplot(2,1,2);
plot(frequency_bpsk,magnitude_bpsk(1:N/2)*2/N);
axis([1700,2300,0,0.8]);
title('扩频后调制信号频谱') ;
xlabel('Hz');


%%
%5.仿真经过awgn信号传输后，扩频后信号时域及频域的变化
%awgn信道模拟了真实的信道，为传输信号增加了高斯白噪声

m_modulation_bitlist_bpsk_awgn=awgn(m_modulation_bitlist_bpsk, 3,'measured');
y_bpsk_awgn=fft(m_modulation_bitlist_bpsk_awgn, N);
magnutide_bpsk_awgn=abs(y_bpsk_awgn);
frequency_bpsk=(1:N/2)*100000/N;


%扩频调制信号经信道后频谱分析
figure(5)
subplot(2,1,1)
plot(ts,m_modulation_bitlist_bpsk,ts,m_modulation_bitlist_bpsk_awgn);
axis([0, 0.005, -1.2, 1.2]);
title('扩频后经加噪过信道后的信号时域波形');

subplot(2,1,2)
plot(frequency_bpsk,magnutide_bpsk_awgn(1:N/2)*2/N);
title('扩频后经信道调制信号频谱');
axis([1700, 2300, 0, 0.8]);
xlabel('Hz');


%%
%6.对比经信道前后两种信号的频谱变化
figure(6)

subplot(2,1,1)
plot(frequency_bpsk, magnitude_bpsk(1:N/2)*2/N);
axis([0,4000,0,0.04]);
title('扩频后调制信号频谱');
xlabel('Hz');

subplot(2,1,2)
plot(frequency_bpsk, magnutide_bpsk_awgn(1:N/2)*2/N);
axis([0,4000,0,0.04]);
title('扩频后施加高斯噪声调制信号频谱');
xlabel('Hz');

%%
%7.不加窄带强干扰
%接收机与本地恢复载波相乘，观察仿真时域波形,观察其频谱变化

bitlist_recovery=m_modulation_bitlist_bpsk_awgn.*cos(2*pi*fs*ts);
frequency_recovery=(1:N/2)*100000/N;
y_bitlist_recovery_fft=fft(bitlist_recovery,N);
magnitude_recovery=abs(y_bitlist_recovery_fft);


figure(7)
subplot(2,1,1)
plot(ts,bitlist_recovery);
axis([0.03, 0.055, -1.5, 1.5]);
title('扩频后接受信号乘以恢复载波')
xlabel('t');

subplot(2,1,2)
plot(frequency_recovery, magnitude_recovery(1:N/2)*2/N);
axis([0, 5000, 0, 0.5]);
xlabel('Hz');
title('扩频后乘以恢复载波后信号频谱');


%%
%8.仿真观察信号经凯萨尔窗低通滤波后的频谱


fp=100;
fc=200;
as=100;
ap=1;%衰减1db
fsw=22000;
wp=2*fp/fsw;
wc=2*fc/fsw;
%求低通滤波器阶数
Nw=ceil((as-7.95)/(14.36*(wc-wp)/2))+1;
beta=0.1102*(as-8.7);
window=kaiser(Nw+1,beta);
%返回截止频率为wc的Nw阶低通滤波器系数向量
b=fir1(Nw,wc,window);
%求频率响应
bs=abs(freqz(b,1,N,fsw))';
magnitude_low_frequency_re=bs.*magnitude_recovery;
%观察经过低通滤波器后无扩频与扩频系统的时域波形
%对扩频系统做iift变换
y_bitlist_recovery_ifft=real(ifft(bs.*y_bitlist_recovery_fft,N));


figure(8)
subplot(2,1,1)
plot(frequency_recovery,magnitude_low_frequency_re(1:N/2)*2/N);
axis([0,200/7,0,0.2]);
title('扩频信号经过低通滤波后频域波形');
xlabel('Hz');

tm=(1:N)/N*4;
subplot(2,1,2)
plot(tm,y_bitlist_recovery_ifft);
xlabel('t');
title('扩频后经过低通滤波后时域波形');

%%
%9，对扩频系统进行解扩，观察其时域频域

%扩频信号乘以解扩序列
bitlist_despread=rectpulse(bipolar_m_sequ,1000);
y_bitlist_recovery_relieve=bitlist_despread.*y_bitlist_recovery_ifft(1:750000);
y_bitlist_recovery_relieve_fft=fft(y_bitlist_recovery_relieve,N);
magnitude_bitlist_recovery_relieve=abs(y_bitlist_recovery_relieve_fft);


figure(9)
subplot(2,1,1)
plot(ts(1:750000),y_bitlist_recovery_relieve(1:750000));
xlabel('t');
axis([0,4,-0.5,0.5]);
title('解扩后信号波形');

subplot(2,1,2)
plot(frequency_recovery,magnitude_bitlist_recovery_relieve(1:N/2)*2/N);
axis([0,500,0,0.2]);
title('解扩后信号频谱');
xlabel('Hz');

%%
%11对解扩信号进行采样，判决

for i=1:750
    K=i*1000-500;
    recovery_bitlist(i)=y_bitlist_recovery_relieve(K);
end
figure(10)
subplot(2,1,1)
stem(recovery_bitlist);
title('解扩信号采样');
for i=1:750
    if recovery_bitlist(i)>0.2
        recovery_bitlist(i)=1;
    elseif recovery_bitlist(i)<-0.2
        recovery_bitlist(i)=-1;
    else
        recovery_bitlist(i)=0;
    end
end
for i=1:50
    K=15*i-14;
    if recovery_bitlist(K)==0
        recovery_bitlist(K)=recovery_bitlist(K+4);
    end
end
for i=1:348
    if recovery_bitlist(i)==0
        recovery_bitlist(i)=recovery_bitlist(i+2);
    end
end
for i=1:50
    output_bitlist(i)=recovery_bitlist(i*15-6);
    if output_bitlist(i)==0
        output_bitlist(i)=output_bitlist(i)+1;
    end
    output_bitlist(i)=(1-output_bitlist(i))/2;
end



subplot(2,1,2)
stem(output_bitlist);
axis([0,50,-0.5,1.5])
title('判决后的最终信号');


figure(11)
subplot(2,1,1)
t=1:50;
stem(t,random_bitlist);
axis([0,50,-0.5,1.5])
title('扩频前二进制比特流');

subplot(2,1,2)
stem(random_bitlist);
axis([0,50,-0.5,1.5])
title('解扩信号采样');

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
function[mseq]=m_sequence(coefficient)
n=length(coefficient);
N=2^n-1; 
register=[zeros(1,n-1) 1]; %移位寄存器的初始状态
mseq(1)=register(n); %m序列的第一个输出码元
for i=2:N
newregister(1)=mod(sum(coefficient.*register),2);
for j=2:n
newregister(j)=register(j-1); 
end;
register=newregister; 
mseq(i)=register(n);
end