%% 
% Sound OFDM Program
% 2011/11/23
% Tom Wada
% 51 subcarrier SOUND OFDM PROGRAM
%

clear;

fprintf('\n');
fprintf('********   SOUND  OFDM  ******** \n');
Fs=44100; % Sampling Frequency=44.1KHz
f0=Fs/1024; % 1KFFT
upconv=40*f0; 
% 40=1.72KHz, 110=4.74KHz, 180=7.75KHz, 250=10.77KHz
centerfreq=upconv; % Center Frequency
fprintf('** Center Frequency = %6.2f KHz ** \n',centerfreq/1000);
%% PILOT GENERATION
rand('seed',0); % Seed for fixed random data
Pdata=floor(rand(51,1)*2); %Pilot Data 
pilot=sqrt(4/3)*(1-2*(Pdata)); %BPSK modelation for Pilot、and multiply squrt(4/3) for boost

%% Sending Data Generation
data0=floor(rand(51,30)*2);
data1=floor(rand(51,30)*2); % Random data generation
data=data1*2+data0;
modqpsk=1/sqrt(2)* [1+1i, -1+1i, 1-1i, -1-1i];
d=modqpsk(1+data); %データをBPSK変調

%% Pilot and Data are combined
dwithP=[pilot,d(:,1:3),pilot,d(:,4:6),pilot,d(:,7:9),pilot,d(:,10:12),pilot,d(:,13:15),pilot,d(:,16:18),pilot,d(:,19:21),pilot,d(:,22:24),pilot,d(:,25:27),pilot,d(:,28:30)];

%% Map to OFDM subcarriers
X=zeros(1024,40); % for 1KFFT and 40 OFDM symbols
X(1:26,:)=dwithP(26:51,:);
X(1000:1024,:)=dwithP(1:25,:);

%constellation
figure(1);
subplot(2,2,1);
plot(X([1:26 1000:1024],:),'.');
axis([-2 2 -2 2]);
title('Constellation at Sender');

%% OFDM modeulation by IFFT
y=128*ifft(X);
%Adding Guard Interval of 512 points
yg=[y(513:1024,:);y];
% 2 dimension data to 1 dimension data
yg=reshape(yg,(1024+512)*40,1);

%% UP CONVERSION
n=1:61440;
LO_cos=(cos(2*pi*centerfreq*n/Fs))';
LO_sin=(-sin(2*pi*centerfreq*n/Fs))';
yy=real(yg).*LO_cos + imag(yg).*LO_sin;

yy=yy*0.3;% amplitude adjustment  should be inside of [-1 to 1]

subplot(2,2,2);
plot(yy)
axis([1 61440 -2 2]);
title('Sounde OFDM waveform');

%% DEMO
sound(yy,Fs);

%% NOISE ADDITION and Multipath
noise =randn((1024+512)*40,1)*0.2;
%Noise
rv=(yy+noise);
%Multipath
delay=100;
rv=yy+0.1*[zeros(delay,1);yy(1:(1024+512)*40-delay,1)]*exp(1i*pi/4)+noise;

rvwv=rv;
%% Spectrum Analysis for received OFDM symbol 
 
Y=fft(rvwv,10000);
Powerspectol = Y .* conj(Y) / 10000;
f = Fs / 10000 / 1000 * (0:5000);
subplot(2,2,3);
plot(f, Powerspectol(1:5001));
title('Power spectral density')
xlabel('Frequency (KHz)')

%% LOW PASS FILTER DESIGN FOR DOWN CONVERSION
Fs = 44100;  % Sampling Frequency

Fpass = 1100;            % Passband Frequency
Fstop = 2200;            % Stopband Frequency
Dpass = 0.057501127785;  % Passband Ripple
Dstop = 0.031622776602;  % Stopband Attenuation
dens  = 20;              % Density Factor

% FIRPMORD function
[N, Fo, Ao, W] = firpmord([Fpass, Fstop]/(Fs/2), [1 0], [Dpass, Dstop]);
% b is coefficient of the LPF
b  = firpm(N, Fo, Ao, W, {dens});

%% DOWN CONVERSION
n=1:61440;
LO_cos=(cos(2*pi*centerfreq*n/Fs))';
LO_sin=(-sin(2*pi*centerfreq*n/Fs))';
% apply LPF
rr = filter(b,1,rvwv.*LO_cos) + 1i*filter(b,1,rvwv.*LO_sin) ;

%% Received signal is reshaped into 2 dimension and DEMODULATION by FFT
ra=reshape(rr,1024+512,40);
ra=ra(513:512+1024,:);
XX=1/512*fft(ra);

%
dd=XX([1000:1024 1:26],:); % DATA and PILOT extraction
da=XX([1000:1024 1:26],1:4:37); %ONLY PILOT extraction
% Channel estimation
hh=[da(:,1)./pilot,da(:,2)./pilot,da(:,3)./pilot,da(:,4)./pilot,da(:,5)./pilot,da(:,6)./pilot,da(:,7)./pilot,da(:,8)./pilot,da(:,9)./pilot,da(:,10)./pilot];
%
%% EQUALIZATION
%パイロットシンボルは1,5,9,13,17,21,25,29,33,37であり
%パイロットシンボルから推定したチャネル伝達関数CTFを用いて、引き続く3シンボルデータをイコライズする。
deq=[];
for n=1:40
    deq=[deq, dd(:,n)./hh(:,ceil(n/4)),];
end
%
subplot(2,2,4);
plot(deq,'.');
axis([-2 2 -2 2]);
title('Constellation after EQUZ|ALIZE');

d0=deq(:,[2:4,6:8,10:12,14:16,18:20,22:24,26:28,30:32,34:36,38:40]);
%% RECOVER DATA
%rdata= (-0.5*sign(real(d0))+0.5)+(-1*sign(imag(d0))+1);
rdata0=-0.5*sign(real(d0))+0.5;
rdata1=-0.5*sign(imag(d0))+0.5;

%% COMPARE SEND DATA and RECEIVED DATA
TotalData=51*30*2;
diff0=sum(sum(abs(sign(data0-rdata0))));
diff1=sum(sum(abs(sign(data1-rdata1))));
fprintf('*** Bit Error Rate *** = %6.4f\n',(diff0+diff1)/TotalData);