lms的频域滤波函数如下：
function [h,y,e] = flmsFunc(x,d,u,N)
%****************** variables *************************
% x : input data
% d : input reference data
% u input step
% N input length of filter
% h output filter coefficient% y output filter data
% e output error
%******************************************************
M = length(x);
y = zeros(1,M);
h = zeros(1,N);
e = zeros(1,M);
for k = 1:fix(M/N) - 1
procx = fft(x((k - 1)*N + 1: (k + 1)*N));
H = fft([h,zeros(1,N)]);
prcox1 = real(ifft(procx.*H));
y(k*N + 1: (k + 1)*N) = prcox1(N + 1:2*N);
e(k*N + 1: (k + 1)*N) = d(k*N + 1: (k + 1)*N) - y(k*N + 1: (k + 1)*N);
E = fft([zeros(1,N),e(k*N + 1: (k + 1)*N)]);
prcox2 = real(ifft(E.*conj(procx)));
v = prcox2(1:N);
h = h + 2*u*v;
end 高斯信道加噪函数：
function [iout,qout] = comb (idata,qdata,attn)
%****************** variables *************************
% idata : input Ich data
% qdata : input Qch data
% iout output Ich data
% qout output Qch data
% attn : attenuation level caused by Eb/No or C/N
%******************************************************
iout = randn(1,length(idata)).*attn;
qout = randn(1,length(qdata)).*attn;
iout = iout+idata(1:length(idata));
qout = qout+qdata(1:length(qdata));
匹配滤波器：
function [iout, qout] = compconv(idata, qdata, filter)
% ****************************************************************
%idata: ich data sequcence
%qdata: qch data sequcence
%filter : filter tap coefficience
% ****************************************************************
iout = conv(idata,filter);
qout = conv(qdata,filter);
生成高斯滤波系数：
function [xh] = gaussf(B,irfn,ipoint,sr,ncc)
%****************************************************************
%irfn: Number of symbols to use filtering
%ipoint: Number of samples in one symbol
%sr: symbol rate
% B: filter coeficiense
%ncc;: 1 -- transmitting filter 0 -- receiving filter
%****************************************************************
tr = sr ;
n = ipoint .* irfn;
mid = ( n ./ 2 ) + 1;
fo=B/sqrt(2*log(2));
for i = 1 : n
icon = i - mid;
ym = icon;
xt=1/2*(erf(-sqrt(2/log(2))*pi*B*(ym/ipoint-1/2)/tr)+erf(sqrt(2/log(2))*pi*B*(ym/ipoint+1/2)/tr));
if ncc == 0 % in the caseof receiver
xh( i ) = xt ;
elseif ncc == 1 % in the caseof transmitter
xh( i ) = xt;
else
error('ncc error');
end
end 下采样函数：
function [out] = oversamp( indata, nsymb , sample)
%****************** variables *************************
% indata : input sequence
% nsymb : Number of symbols
% sample : Number of oversample
% *****************************************************
out=zeros(1,nsymb*sample);
out(1:sample:1+sample*(nsymb-1))=indata;
主函数，包含GMSK调制解调以及flms滤波过程：
function gmskModDemodFlms()
berMatrix = [];
ebn0Matrix = -10:2:10;
sr=256000.0; % Symbol rate
ml=1; % ml:Number of modulation levels
br=sr.*ml; % Bit rate
nd = 1000; % Number of symbols that simulates in each loop
IPOINT=8; % Number of oversamples
%************************* Filter initialization ***************************
irfn=21; % Number of taps
B=0.15*sr;
B2=0.4*sr;
[xh] = gaussf(B,irfn,IPOINT,sr,1); %Transmitter filter coefficients
[xh2] =gaussf(B2,irfn,IPOINT,sr,0); %Receiver filter coefficients
u = 0.0001;
N = 5;
for ebn0 = ebn0Matrix%Eb/N0
%******************** START CALCULATION *************************************
nloop=100; % Number of simulation loops
noe = 0; % Number of error data
nod = 0; % Number of transmitted data
for k=1:nloop
%*************************** Data generation ********************************
data1=rand(1,nd.*ml)>0.5; % rand: built in function
%*************************** GMSK Modulation ********************************
data11=2*data1-1;
data2=oversamp(data11,length(data11),IPOINT);
data3=conv(data2,xh); % NEW forGMSK
th=zeros(1,length(data3)+1);
ich2=zeros(1,length(data3)+1);
qch2=zeros(1,length(data3)+1);
for ii=2:length(data3)+1
th(1,ii)=th(1,ii-1)+pi/2*data3(1,ii-1)./IPOINT;
end
ich2=cos(th);
qch2=sin(th);
%************************** Attenuation Calculation ***********************
spow=sum(ich2.*ich2+qch2.*qch2)/nd; % sum: built in function
attn=0.5*spow*sr/br*10.^(-ebn0/10);
attn=sqrt(attn); % sqrt: built in function
%********************** Fading channel **********************
%[ifade,qfade]=sefade2(data2,qdata1,itau,dlvl1,th1,n0,itnd1,now1,length(data2),fftlen2,fstp,fd,flat);
%********************* Add White Gaussian Noise (AWGN) **********************
[ich3,qch3]= comb(ich2,qch2,attn);% add white gaussian noise
[ich4,qch4] = compconv(ich3,qch3,xh2);
syncpoint =irfn*IPOINT-IPOINT/2+1;
ich5=ich4(syncpoint:IPOINT:length(ich4));
qch5=qch4(syncpoint:IPOINT:length(qch4));
d1 = ich5 - randn(size(ich5))*attn*0.001;
d2 = qch5 - randn(size(qch5))*attn*0.001;
[~,ich5,~] = flmsFunc(ich5,d1,u,N);
[~,qch5,~] = flmsFunc(qch5,d2,u,N);
%**************************** GMSK Demodulation *****************************
demoddata2(1,1)=-1;
for kk=3:2:nd*ml+1
demoddata2(1,kk)=ich5(1,kk)*qch5(1,kk-1)*cos(pi*(kk))>0;
end
for n=2:2:nd*ml+1
demoddata2(1,n)=ich5(1,n-1)*qch5(1,n)*cos(pi*(n))>0;
end
[demodata]=demoddata2(1,2:nd*ml+1);
%************************** Bit Error Rate (BER) ****************************
noe2=sum(abs(data1-demodata)); % sum: built in function
nod2=length(data1); % length: built in function
noe=noe+noe2;
nod=nod+nod2;
end
%********************** result ***************************
ber = noe/nod;
berMatrix = cat(1,berMatrix,ber);
end