run.asv

很好的ofdm文档

B= 4; % # of bits
C= 16; % # of points in constellation
S= 256; % # samples of handel
L=16; % length of channel
noise=.0001; % noise amplitude
sN=0; % # TD spikes
N=L; % length of cyclic pad
F=128; % length of blocks (1/2 # of subchannels)
M=constbuilder(C); % constellation matrix

s=abs(randn(L,1));  % random channel
s=s./sum(s);

load handel
y=y(1:S);
x=a2d(y,B);
z=constmap(x,M);
a=series2parallel(z,F);
b=mirror(a);
c=ifft(b); 
d=cyclicpad(c,N);
e=matlaberr(parallel2series(d)); 

aN = randn(length(conv(e,s)),1)*noise;

f=channel(e,s,aN,sN);
f=f(1:length(e));
g=series2parallel(f,2*F+N);
h=decyclicpad(g,N);%
k=fft(h); 
m=filterchannel(k,s);
n=demirror(m);
n(1,:)=a(1,:);
o=parallel2series(n); 
p=approximate(o,M);
q=deconstmap(p,M);
r=d2a(q,B);

% figure(1)
% stem(y/max(abs(y)));title('input speech');
% axis([1 S -1 1]);
% grid on
%     
% figure(2)
% stem(d2a(x,B));ylabel('quantized input');
% axis([1 S -1 1]);
% grid on
% 
% figure(3)
% plot(real(M(:,length(M(1,:)))),imag(M(:,length(M(1,:)))),'.b');
% grid on
% xlabel('REAL');
% ylabel('IMAG');
% axis([-1.5 1.5 -1.5 1.5]);
% 
% figure(4)
% Q=length(d(:,1));
% stem(1:Q-N,matlaberr(d(1:Q-N,1)),'b');ylabel('cylic padded TD input');
% hold
% stem(1:N,matlaberr(d(1:N,1)),'r');
% stem(Q-N+1:Q,matlaberr(d((Q-N+1):Q,1)),'r');
% hold
% 
% 
% figure(5)
% %stem(e);ylabel('series TD input');
% Q=2*(2*F+N);
% stem(1:Q/2,e(1:Q/2));
% hold
% stem(Q/2+1:Q,e(Q/2+1:Q),'g');
% stem(Q/2-L/2+1/2:Q/2+L/2,s,'r');
% hold
% 
% figure(7)
% plot(real(o),imag(o),'*r');
% hold;
% plot(real(M(:,length(M(1,:)))),imag(M(:,length(M(1,:)))),'.b');
% ylabel('Received Constellation Points');
% hold;
% grid;
% axis([-1.5 1.5 -1.5 1.5]);

% figure(6)
% stem(fftshift(abs(fft(s,256))));ylabel('Transfer Function of Channel');
% axis([0 2*F 0 1]);
% 
% figure(8)
% stem(f);ylabel('TD output of channel');
% 
% figure(9)
% stem(h);ylabel('cyclic pad removed');
% 
% figure(10)
% stem(fftshift(abs(1./fft(s,256))));ylabel('transfer function of equalizer');
%   
% figure(11)
% stem(r);ylabel('received/quantized speech'); 
% 
% figure(12)
% stem(1:128,abs(b(1:length(b)/2,1)),'m')
% hold
% stem(129:256,abs(b(length(b)/2+1:length(b),1)),'k')
% hold

ps=(norm(e))^2;
pn=(norm(aN))^2;
snr=10*log10(ps/pn);
biterr=sum(abs(x-q));

function y = constbuilder(N)
b=ceil(log2(N));
p=zeros(N,b);
for k=0:(N-1)
        u=dec2bin(k);
        w=zeros(1,length(u));
        for i=1:length(u)
            w(i)=str2num(u(i));
        end
        v=[zeros((b-length(w)),1)',w];
        p(k+1,:)=v;
end
n=0:N-1;
    if N<=16
        a=exp(j*2*pi/N.*n');
        M=[p,a];
    else
        a=exp(j*2*pi/16.*(0:15)');
        if N<=24
            b=0.5*exp(j*2*pi/(N-17).*(0:N-17)');
            M=[p,[a;b]];
        else
            b=0.66*exp(j*2*pi/8.*(0:8)'+j*pi/16);
            if N<=32
                c=0.33*exp(j*2*pi/(N-26).*(0:N-26)');
                M=[p,[a;b;c]];
            else
                a=exp(j*2*pi/N.*n');
                M=[p,a];
            end
        end
    end
if N==16
    a=[.354+.354*j;sqrt(2)/2;sqrt(2)/2*j;-.354+.354*j;-sqrt(2)/2*j;.354-.354*j;-.354-.354*j;-sqrt(2)/2];
    b=exp(j*2*pi/8.*(0:7)');
    M=[p,[a;b]];
else
end    
y=matlaberr(M);



function y = a2d(x,b)
% a2d(x) - converts discrete time, continuous amplitude signal to
%          bitstream (b bits per sample)
%  x - original signal
%  b - number bits per sample 

N = length(x);
qinterval = 2*max(abs(x))/((2^b)-1);
if b<=1;
    for i=1:N
        if x(i)<=0
            bits(i)=0;
        else
            bits(i)=1;
        end
    end
else
    q = quant(x-min(x),qinterval);
    bits=zeros(N*b,1);
    
    for k=0:((2^b)-1)
        u=dec2bin(k);
        for i=1:length(u)
            w(i)=str2num(u(i));
        end
        v=[zeros((b-length(w)),1)',w];
        p(k+1,:)=v;
    end
    for i=0:N-1
        for k=0:((2^b)-1)
            if q(i+1)==k*qinterval
                bits(b*i+1:b*i+b)=p(k+1,:);
                break
            else
            end
        end
    end
end
y=bits;


function y = constmap(x,M)
% constmap(x) - maps constellation onto bitstream (2 bits/symbol)
%   x - bitstream

m=length(M(1,:))-1;
N = length(x);
newN=m*round(N/m);
x=x(1:newN);
for i=1:m
    newvector(:,i)=x(i:m:newN-(m-i));
end

alpha=zeros(newN/m,1);
for i=1:newN/m
    for k=1:length(M(:,1))
        if (newvector(i,:)==M(k,1:m))
        alpha(i)=M(k,m+1);
        break
    else
    end
end
end
y=alpha;

function y = series2parallel(x,N)

L=length(x);
q=floor(L/N);
newvec=zeros(N,q);
for i=1:q
    newvec(1:N,i)=x((1+(i-1)*N):i*N);
end
y=newvec;


function y = mirror(X)
R=length(X(1,:));
L=length(X(:,1));
if isreal(X(1,:))
    X=X;
else
    X(1,:)=zeros(1,R);
end
    newvec=zeros(2*L,R);
    newvec(1:L,:)=X;
    newvec(L+1,:)=X(1,:);
    newvec(L+2:2*L,:)=conj(flipud(X(2:L,:)));
y=newvec;


function y=cyclicpad(X,L)
N=length(X(:,1));
Y=[X(N-L+1:N,:);X];
y=Y;


function y = matlaberr(X)
R=length(X(1,:));
L=length(X(:,1));
for m=1:R
    for n=1:L
        if abs(real(X(n,m)))<.000001
            a=0;
        else
            a=real(X(n,m));
        end
        if abs(imag(X(n,m)))<.000001
            b=0;
        else 
            b=imag(X(n,m));
        end
        X(n,m)=a+j*b;
    end
end
y=X;


function y = parallel2series(X)
L=length(X(:,1));
R=length(X(1,:));
newvec=zeros(L*R,1);
for i=0:(R-1)
    newvec(1+i*L:(1+i)*L,1)=X(:,i+1);
end
y=newvec;

function y=channel(x,h,aN,sN)

a=conv(x,h);
N=length(a);
w=a+aN;
m=max(abs(x));
for k=1:sN
    w(round(abs(randn(1,1))*(N-1))+1)=m; 
end
y=w;

function y=decyclicpad(X,L)
N=length(X(:,1));
Y=X(L+1:N,:);
y=Y;


function y=filterchannel(X,h)
N=length(X(:,1));
W=length(X(1,:));

H=fft(h,N);

invH=1./H;
for k=1:W
    y(:,k)=invH.*X(:,k);
end
y;


function y = demirror(X)
R=length(X(1,:));
L=length(X(:,1));
if rem(L,2)==1
    y=[ones(1,R);X(1:L/2,:)];
else
    y=X(1:L/2,:);
end

function y = approximate(x,M)
% approximate(x,M) - approximates noisy symbol vector x with valid symbols
% in M
%   x - noisy symbol vector
%   M - matrix representing constellation mapping

%M=[0 0 1;0 1 j;1 0 -1;1 1 -j];
W=length(M(1,:));
L=length(M(:,1));
N = length(x);
a = M(:,W);
alpha=zeros(N,1);

for i=1:N
    mindist=100;
    symnum=0;
    for k=1:L
        if abs(a(k)-x(i))<=abs(mindist)
            mindist=abs(a(k)-x(i));
            alpha(i)=a(k);
        else
        end
    end
end        
y=alpha;

function y = deconstmap(x,M)
% deconstmap(x) - maps constellation into bitstream (2 bits/symbol)
%   x - bitstream

%M=[0 0 1;0 1 j;1 0 -1;1 1 -j];
m=length(M(1,:))-1;   % # of bits/symbol
N = length(x);
bits=zeros(N*m,1);
for i=0:N-1
    for k=1:length(M(:,1))
        if (x(i+1)==M(k,m+1))
        bits(1+(m*i):m*(i+1))=M(k,1:m);
        break
    else
    end
end
end

y=bits;

function y = deconstmap(x,M)
% deconstmap(x) - maps constellation into bitstream (2 bits/symbol)
%   x - bitstream

%M=[0 0 1;0 1 j;1 0 -1;1 1 -j];
m=length(M(1,:))-1;   % # of bits/symbol
N = length(x);
bits=zeros(N*m,1);
for i=0:N-1
    for k=1:length(M(:,1))
        if (x(i+1)==M(k,m+1))
        bits(1+(m*i):m*(i+1))=M(k,1:m);
        break
    else
    end
end
end

y=bits;