Matlab Code

%% 

f = 0:1/1024:(152-1)/1024;

N = 41;

n = 0:1:N-1;

C=41;

c=0:1:C-1;

%x =- 2*sin(2*pi*2*n)-0.3*sin(2*pi*0.2*n);%%%%%discretized low frequency function

x=-.25*n+2+.99*n;

x2=(-1).^c;%%low to high frequency conversion function

t=x.*x2;%%High frequency function response

figure

subplot(3,1,1),stem(x),title('señal original muestreada')

%subplot(3,1,2),stem(x2),title('señal muestreada en magnitud 1 a -1')

%subplot(3,1,3),stem(t),title('señal convertida a altas frecuencias muestreada ')

%freqz(t)

X=fft(x,1024)

X2=fft(x2,1024)

T=fft(t,1024)

figure

subplot(3,1,1),plot(f,abs(X(1:152))),title('señal original muestreada en dominio de frecuencia')

subplot(3,1,2),plot(f,abs(X2(1:152))),title('señal muestreada en dominio de frecuencia de magnitud 1 a -1')

subplot(3,1,3),plot(f,abs(T(1:152))),title('señal convertida a altas frecuencias, muestreada en dominio de frecuencia')

%% caso 1

%%

%% 

t = 0:0.01:2;

%h1=firgr(80,[0,wp,ws,1]   ,   [1,1,0,0 ]);

%x1=downsample(x,2);

mt = cos(2*pi*t)

x = cos((2*pi) * (20*t)) .* mt

figure

plot(x)

%%

I = imread('cameraman.tif');

I = im2double(I); %Convierte la imagen a valores reales

T = dctmtx(8); %Divide la imagen en bloques 8x8

dct = @(block_struct) T * block_struct.data * T';

B = blockproc(I,[8 8],dct);

mask = [1   1   1   1   0   0   0   0

               1   1   1   0   0   0   0   0

               1   1   0   0   0   0   0   0

               1   0   0   0   0   0   0   0

               0   0   0   0   0   0   0   0

               0   0   0   0   0   0   0   0

               0   0   0   0   0   0   0   0

               0   0   0   0   0   0   0   0];

B2 = blockproc(B,[8 8],@(block_struct) mask .* block_struct.data);

invdct = @(block_struct) T' * block_struct.data * T;

I2 = blockproc(B2,[8 8],invdct);

imshow(I), figure, imshow(I2)

%% 

%Caso 2

t = 0:0.01:2;

n= 0:0.01:2;

mt = cos(2*pi*0.5*t);

Mt = cos(2*pi*0.5*n/8000);

figure

x = cos((2*pi) * (22/3*t)) .* mt.*cos(2*pi*2*t);

X = cos((2*pi) * (22/3*n/8000)) .* Mt.*cos(2*pi*2*n/4);

%%m(t)

plot(x),title('x(t)')

%plot(X),title('x(nT)')

%Figura compuesta

plot(t,x,'g',t,X,'r');

%%

%Caso 1

t = 0:0.01:2;

f = 0:1/1024:(152-1)/1024;

n =  0:0.01:2;

mt = cos(2*pi*t);

Mt = cos(2*pi*n/(8000));

x = cos((2*pi) * (18000/4*t)) .* mt;

X = cos((2*pi)*(18000/4*n/(8000))).* Mt;

%%m(t)

%m = cos((2*pi) * t);

%Figura compuesta

figure

stem(X)

XX=fft(x,1024);

XXX=fft(X,1024);

subplot(2,1,1),plot(f,abs(XX(1:152)))

subplot(2,1,2),plot(f,abs(XXX(1:152)))

%%

%Caso 2

t = 0:0.01:2;

n= 0:0.01:2;

mt = cos(2*pi*0.5*t);

Mt = cos(2*pi*0.5*n/4);

figure

x = cos((2*pi) * (22*t)) .* mt;

X = cos((2*pi) * (22*t)) .* Mt;

%%m(t)

plot(x),title('x(t)')

%plot(X),title('x(nT)')

%Figura compuesta

%plot(t,x,'g',t,X,'r');

%% Prubas de tesis de maestria Nicthe Mtz 26/Oct/2015

%-------------------------------------------------------------------------------

%------------------------------------------------------------------------------

r=0.05;

t=41;

T= 0:1:T-1;

Ten=(1-r)/2T:(1+r)/2T;

GRC=Ten/2(1+cos(pi*Ten/r(f-(1-r)/2Ten)));

G=fft(GRC,1024);

%% 

% test_orthogonality.m

% To check the orthogonality among some sinusoidal signals

% with different frequencies/phases

clear, clf

T=1.6; ND=1000; nn=0:ND; ts=0.002; tt=nn*ts; % Time interval

Ts = 0.1; M = round(Ts/ts); % Sampling period in continuous/discrete-tim

nns = [1:M:ND+1]; tts = (nns-1)*ts; % Sampling indices and times

ks = [1:4 3.9 4]; tds = [0 0.1 0.9 0 0.99 0]; % Frequencies and delays

K = length(ks);

for i=1:K

k=ks(i); td=tds(i); x(i,:) = exp(j*2*pi*k*(tt-td)/T);

if i==K, x(K,:) = [x(K,[302:end]) x(K-3,[1:301])]; end

subplot(K,2,2*i-1), plot(tt,real(x(i,:))),

hold on, plot(tt([1 end]),[0 0],'k'), stem(tts,real(x(i,nns)),'.')

end

N = round(1.6/0.1); 

Xk = fft(.').'; kk = 0:N-1;

for i=1:K,

k=ks(i); td=tds(i); subplot(K,2,2*i), stem(kk,abs(Xk(i,:)),'.');

end

%% 

%% 

% OFDM_basic.m

clear all

NgType=1; % NgType=1/2 for cyclic prefix/zero padding

if NgType==1, nt='CP'; elseif NgType==2, nt='ZP'; end

Ch=0; % Ch=0/1 for AWGN/multipath channel

if Ch==0, chType='AWGN'; Target_neb=100; else chType='CH'; Target_neb=500; end

figure(Ch+1), clf

PowerdB=[0 -

8 -17 -21 -25]; % Channel tap power profile dB%%

Delay=[0 3 5 6 8]; % Channel delay sample

Power=10.^(PowerdB/10); % Channel tap power profile linear scale

Ntap=length(PowerdB); % Chanel tap number

Lch=Delay(end)+1; % Channel length-----------------------------------------------------

Nbps=4; M=2^Nbps; % Modulation order=2/4/6 for QPSK/16QAM/64QAM

Nfft=64; % FFT size

Ng=Nfft/4; % GI (Guard Interval) length (Ng=0 for no GI)

Nsym=Nfft+Ng; % Symbol duration

Nvc=Nfft/4; % Nvc=0: no VC (virtual carrier)

Nused=Nfft-Nvc;

EbN0=[0:5:30]; % EbN0

N_iter=1e5; % Number of iterations for each EbN0

Nframe=3; % Number of symbols per frame

sigPow=0; % Signal power initialization

file_name=['OFDM_BER_' chType '_' nt '_' 'GL' num2str(Ng) '.dat'];

fid=fopen(file_name, 'w+');

norms=[1 sqrt(2) 0 sqrt(10) 0 sqrt(42)]; % BPSK 4-QAM 16-QAM

for i=0:length(EbN0)

randn('state',0); rand('state',0); Ber2=ber(); % BER initialization

Neb=0; Ntb=0; % Initialize the number of error/total bits

for m=1:N_iter

% Tx______________________________________________________________

X= randint(1,Nused*Nframe,M); % bit: integer vector

Xmod= qammod(X,M,0,'gray')/norms(Nbps);

if NgType=2, x_GI=zeros(1,Nframe*Nsym);

elseif NgType==2, x_GI= zeros(1,Nframe*Nsym+Ng);

% Extend an OFDM symbol by Ng zeros

end

kk1=[1:Nused/2]; kk2=[Nused/2+1:Nused]; kk3=1:Nfft; kk4=1:Nsym;

for k=1:Nframe

if Nvc=0, X_shift= [0 Xmod(kk2) zeros(1,Nvc-1) Xmod(kk1)];

else X_shift= [Xmod(kk2) Xmod(kk1)];

end

Introduction to OFDM 137

x= ifft(X_shift);--------------------------------------------------------------------

x_GI(kk4)= guard_interval(Ng,Nfft,NgType,x);-------------------------------------------------

kk1=kk1+Nused; kk2= kk2+Nused; kk3=kk3+Nfft; kk4=kk4+Nsym;----------------------------------

end

if Ch==0, y= x_GI; % No channel

else % Multipath fading channel

channel=(randn(1,Ntap)+j*randn(1,Ntap)).*sqrt(Power/2);%----------------------------------------------------

h=zeros(1,Lch); h(Delay+1)=channel; % cir: channel impulse response-------------------------------------------

y = conv(Xn,h);%------------------------------------------------------------------------------------------

end

if i==0 % Only to measure the signal power for adding AWGN noise

y1=y(1:Nframe*Nsym); sigPow = sigPow + y1*y1'; continue;

end

% Add AWGN noise________________________________________________

snr = EbN0(i)+10*log10(Nbps*(Nused/Nfft)); % SNR vs. Eb/N0 by Eq.(4.28)

noise_mag = sqrt((10.^(-snr/10))*sigPow/2);

y_GI = y + noise_mag*(randn(size(y))+j*randn(size(y)));

% Rx_____________________________________________________________

kk1=(NgType==2)*Ng+[1:Nsym]; kk2=1:Nfft;

kk3=1:Nused; kk4=Nused/2+Nvc+1:Nfft; kk5=(Nvc=0)+[1:Nused/2];

if Ch==1

H= fft([h zeros(1,Nfft-Lch)]); % Channel frequency response

H_shift(kk3)= [H(kk4) H(kk5)];

end

for k=1:Nframe

Y(kk2)= fft(remove_GI(Ng,Nsym,NgType,y_GI(kk1)));

Y_shift=[Y(kk4) Y(kk5)];

if Ch==0, Xmod_r(kk3) = Y_shift;

else Xmod_r(kk3)=Y_shift./H_shift; % Equalizer - channel compensation

end

kk1=kk1+Nsym; kk2=kk2+Nfft; kk3=kk3+Nused; kk4=kk4+Nfft;

kk5=kk5+Nfft;

end

X_r=qamdemod(Xmod_r*norms(Nbps),M,0,'gray');

Neb=Neb+sum(sum(de2bi(X_r,Nbps)=de2bi(X,Nbps)));

Ntb=Ntb+Nused*Nframe*Nbps; %[Ber,Neb,Ntb]=ber(bit_Rx,bit,Nbps);

if Neb>Target_neb, break; end

end%

if i==0, sigPow= sigPow/Nsym/Nframe/N_iter;

else

Ber = Neb/Ntb;

fprintf('EbN0=%3d[dB], BER=%4d/%8d =%11.3e\n’, EbN0(i), Neb,Ntb,Ber)

fprintf(fid, '%d\t%11.3e\n’, EbN0(i), Ber);

if Ber<1e-6, break; end

end

end

if (fid =0), fclose(fid); end

plot_ber(file_name,Nbps);

%138 MIMO-OFDM Wireless Communications with MATLAB

Progr am 4.3 Routines for GI (guard interval) insertion , GI removal, and BER p lotting

function y = guard_interval(Ng,Nfft,NgType,ofdmSym)

if NgType==1, y=[ofdmSym(Nfft-Ng+1:Nfft) ofdmSym(1:Nfft)];

elseif NgType==2, y=[zeros(1,Ng) ofdmSym(1:Nfft)];

end

function y=remove_GI(Ng,Lsym,NgType,ofdmSym)

if Ng  =0

if NgType==1, y=ofdmSym(Ng+1:Lsym); % cyclic prefix

elseif NgType==2 % cyclic suffix

y=ofdmSym(1:Lsym-Ng)+[ofdmSym(Lsym-Ng+1:Lsym) zeros (1,Lsym-2*Ng)];

end

else y=ofdmSym;

end

function plot_ber(file_name,Nbps)

EbN0dB=[0:1:30]; M=2^Nbps;

ber_AWGN = ber_QAM(EbN0dB,M,'AWGN');

ber_Rayleigh = ber_QAM(EbN0dB,M,'Rayleigh');

semilogy(EbN0dB,ber_AWGN,'r:'), hold on,

semilogy(EbN0dB,ber_Rayleigh,'r-')

a= load(file_name); semilogy(a(:,1),a(:,2),'b–s'); grid on

legend('AWGN analytic','Rayleigh fading analytic', 'Simulation');

xlabel('EbN0[dB]'), ylabel('BER'); axis([a(1,1) a(end,1) 1e-5 1])

function ber=ber_QAM(EbN0dB,M,AWGN_or_Rayleigh)

% Find analytical BER of M-ary QAM in AWGN or Rayleigh channel

% EbN0dB=EbN0dB: Energy per bit-to-noise power[dB] for AWGN channel

% =rdB : Average SNR(2*sigma Eb/N0)[dB] for Rayleigh channel

% M = Modulation order (Alphabet or Constellation size)

N= length(EbN0dB); sqM= sqrt(M);

a= 2*(1-power(sqM,-1))/log2(sqM); b= 6*log2(sqM)/(M-1);

if nargin<3, AWGN_or_Rayleigh='AWGN'; end

if lower(AWGN_or_Rayleigh(1))=='a'

ber = a*Q(sqrt(b*10.^(EbN0dB/10))); % ber=berawgn(EbN0dB,’QAM’,M) Eq.(4.25)

else % diversity_order=1; ber=berfading(EbN0dB,’QAM’,M,diversity_order)

rn=b*10.^(EbN0dB/10)/2; ber = 0.5*a*(1-sqrt(rn./(rn+1))); % Eq.(4.26)

end

function y=Q(x)

% co-error function: 1/sqrt(2*pi) * int_x^inf exp(-t^2/2) dt. % Eq.(4.27)

y=erfc(x/sqrt(2))/2;