diff --git a/lte/phy/lib/phch/test/pcfich_test_mex.c b/lte/phy/lib/phch/test/pcfich_test_mex.c
new file mode 100644
index 000000000..4fff6cc7e
--- /dev/null
+++ b/lte/phy/lib/phch/test/pcfich_test_mex.c
@@ -0,0 +1,157 @@
+/**
+ *
+ * \section COPYRIGHT
+ *
+ * Copyright 2013-2014 The libLTE Developers. See the
+ * COPYRIGHT file at the top-level directory of this distribution.
+ *
+ * \section LICENSE
+ *
+ * This file is part of the libLTE library.
+ *
+ * libLTE is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, either version 3 of
+ * the License, or (at your option) any later version.
+ *
+ * libLTE is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * A copy of the GNU Lesser General Public License can be found in
+ * the LICENSE file in the top-level directory of this distribution
+ * and at http://www.gnu.org/licenses/.
+ *
+ */
+
+#include <string.h>
+#include "liblte/phy/phy.h"
+#include "liblte/mex/mexutils.h"
+
+/** MEX function to be called from MATLAB to test the channel estimator 
+ */
+
+#define ENBCFG  prhs[0]
+#define INPUT   prhs[1]
+#define NOF_INPUTS 2
+
+void help()
+{
+  mexErrMsgTxt
+    ("[cfi] = liblte_pdcch(enbConfig, rxWaveform)\n\n");
+}
+
+/* the gateway function */
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+  int i; 
+  lte_cell_t cell; 
+  pcfich_t pcfich;
+  chest_dl_t chest; 
+  lte_fft_t fft; 
+  regs_t regs;
+  uint32_t sf_idx; 
+  cf_t *input_fft, *input_signal;
+  
+  if (nrhs != NOF_INPUTS) {
+    help();
+    return;
+  }
+    
+  if (mexutils_read_cell(ENBCFG, &cell)) {
+    help();
+    return;
+  }
+  
+  if (mexutils_read_uint32_struct(ENBCFG, "NSubframe", &sf_idx)) {
+    help();
+    return;
+  }
+
+  if (chest_dl_init(&chest, cell)) {
+    fprintf(stderr, "Error initializing equalizer\n");
+    return;
+  }
+
+  if (lte_fft_init(&fft, cell.cp, cell.nof_prb)) {
+    fprintf(stderr, "Error initializing FFT\n");
+    return;
+  }
+  
+  if (regs_init(&regs, cell)) {
+    mexErrMsgTxt("Error initiating regs\n");
+    return;
+  }
+  
+  if (pcfich_init(&pcfich, &regs, cell)) {
+    fprintf(stderr, "Error creating PBCH object\n");
+    return;
+  }
+      
+  /** Allocate input buffers */
+  if (mexutils_read_cf(INPUT, &input_signal) < 0) {
+    mexErrMsgTxt("Error reading input signal\n");
+    return; 
+  }
+  input_fft = vec_malloc(SF_LEN_RE(cell.nof_prb, cell.cp) * sizeof(cf_t));
+  
+  // Set Channel estimates to 1.0 (ignore fading) 
+  cf_t *ce[MAX_PORTS];
+  for (i=0;i<cell.nof_ports;i++) {
+    ce[i] = vec_malloc(SF_LEN_RE(cell.nof_prb, cell.cp) * sizeof(cf_t));
+  }
+  
+  lte_fft_run_sf(&fft, input_signal, input_fft);
+
+  if (nrhs > NOF_INPUTS) {
+    cf_t *cearray; 
+    mexutils_read_cf(prhs[NOF_INPUTS], &cearray);
+    for (i=0;i<cell.nof_ports;i++) {
+      for (int j=0;j<SF_LEN_RE(cell.nof_prb, cell.cp);j++) {
+        ce[i][j] = *cearray;
+        cearray++;
+      }
+    }
+  } else {
+    chest_dl_estimate(&chest, input_fft, ce, sf_idx);    
+  }
+  float noise_power;
+  if (nrhs > NOF_INPUTS + 1) {
+    noise_power = mxGetScalar(prhs[NOF_INPUTS+1]);
+  } else {
+    noise_power = chest_dl_get_noise_estimate(&chest);
+  }
+    
+    
+  uint32_t cfi, distance; 
+  int n = pcfich_decode(&pcfich, input_fft, ce, noise_power,  sf_idx, &cfi, &distance);
+
+  if (nlhs >= 1) { 
+    if (n < 0) {      
+      plhs[0] = mxCreateDoubleScalar(-1);
+    } else {
+      plhs[0] = mxCreateDoubleScalar(cfi);      
+    }
+  }
+  if (nlhs >= 2) {
+    mexutils_write_cf(pcfich.pcfich_d, &plhs[1], 16, 1);  
+  }
+  if (nlhs >= 3) {
+    mexutils_write_cf(pcfich.pcfich_symbols[0], &plhs[2], 16, 1);  
+  }
+  
+  chest_dl_free(&chest);
+  lte_fft_free(&fft);
+  pcfich_free(&pcfich);
+  regs_free(&regs);
+  
+  for (i=0;i<cell.nof_ports;i++) {
+    free(ce[i]);
+  }
+  free(input_signal);
+  free(input_fft);
+  
+  return;
+}
+
diff --git a/matlab/tests/pbch_bler.m b/matlab/tests/pbch_bler.m
new file mode 100644
index 000000000..1802dbb83
--- /dev/null
+++ b/matlab/tests/pbch_bler.m
@@ -0,0 +1,104 @@
+%clear
+% R.1 10 MHz 1 port
+% R.10 10 MHz 2 ports
+% R.4 1.4 MHz 1 port 
+% R.11-2 5 MHz 2 ports
+rmc = lteRMCDL('R.10');
+
+NofPortsTx=2;
+
+SNR_values_db=linspace(-8,-2,4);
+Nrealizations=200;
+enb = struct('NCellID',0,'NDLRB',50,'CellRefP',NofPortsTx,'CyclicPrefix','Normal','DuplexMode','FDD','NSubframe',0);
+
+griddims = lteResourceGridSize(enb); % Resource grid dimensions
+L = griddims(2);    
+        
+cfg.Seed = 8;                  % Random channel seed
+cfg.NRxAnts = 1;               % 1 receive antenna
+cfg.DelayProfile = 'EPA';      % EVA delay spread
+cfg.DopplerFreq = 5;           % 120Hz Doppler frequency
+cfg.MIMOCorrelation = 'Low';   % Low (no) MIMO correlation
+cfg.InitTime = 0;              % Initialize at time zero
+cfg.NTerms = 16;               % Oscillators used in fading model
+cfg.ModelType = 'GMEDS';       % Rayleigh fading model type
+cfg.InitPhase = 'Random';      % Random initial phases
+cfg.NormalizePathGains = 'On'; % Normalize delay profile power 
+cfg.NormalizeTxAnts = 'On';    % Normalize for transmit antennas
+
+cec.PilotAverage = 'UserDefined';     % Type of pilot averaging
+cec.FreqWindow = 9;                   % Frequency window size
+cec.TimeWindow = 9;                   % Time window size
+cec.InterpType = 'linear';             % 2D interpolation type
+cec.InterpWindow = 'Centered';        % Interpolation window type
+cec.InterpWinSize = 1;                % Interpolation window size
+
+rmc.PDSCH.Modulation = '16QAM';
+[waveform,rgrid,info] = lteRMCDLTool(rmc,[1;0;0;1]);
+
+cfg.SamplingRate = info.SamplingRate;
+
+addpath('../../debug/lte/phy/lib/phch/test')
+
+
+error=zeros(length(SNR_values_db),2);
+for snr_idx=1:length(SNR_values_db)
+    SNRdB = SNR_values_db(snr_idx);             % Desired SNR in dB
+    SNR = 10^(SNRdB/20);    % Linear SNR  
+
+    errorReal = zeros(Nrealizations,2);
+    for i=1:Nrealizations
+
+        rxWaveform = lteFadingChannel(cfg,sum(waveform,2));
+
+        %% Additive Noise
+        N0 = 1/(sqrt(2.0*double(enb.CellRefP)*double(info.Nfft))*SNR);
+
+        % Create additive white Gaussian noise
+        noise = N0*complex(randn(size(rxWaveform)),randn(size(rxWaveform)));   
+
+        rxWaveform = noise + rxWaveform;
+        
+       % rxWaveform = downsampled; 
+        
+        % Number of OFDM symbols in a subframe
+        % OFDM demodulate signal
+        rxgrid = lteOFDMDemodulate(enb, rxWaveform);
+
+        % Perform channel estimation
+        [hest, nest] = lteDLChannelEstimate(enb, cec, rxgrid(:,1:L,:));
+      
+        pbchIndices = ltePBCHIndices(enb);
+        [pbchRx, pbchHest] = lteExtractResources( ...
+            pbchIndices, rxgrid(:,1:L,:), hest(:,1:L,:,:));
+
+        % Decode PBCH
+        [bchBits, pbchSymbols, nfmod4, mib, nof_ports] = ltePBCHDecode(enb, pbchRx, pbchHest, nest);
+
+        if (nof_ports ~= NofPortsTx)
+            errorReal(i,1)=1;
+        end
+        
+        [nof_ports2, pbchSymbols2, pbchBits, ce, ce2, pbchRx2, pbchHest2,indices]= liblte_pbch(enb, rxWaveform, hest, nest);
+        if (nof_ports2 ~= NofPortsTx)
+            errorReal(i,2)=1;
+        end
+%         if (errorReal(i,1) ~= errorReal(i,2))
+%             i=1;
+%         end
+    end
+    error(snr_idx,:) = sum(errorReal);
+    fprintf('SNR: %.2f dB\n', SNR_values_db(snr_idx));
+end
+
+if (length(SNR_values_db) > 1)
+    semilogy(SNR_values_db, error/Nrealizations)
+    grid on
+    xlabel('SNR (dB)');
+    ylabel('BLER')
+    legend('Matlab','libLTE')
+    axis([min(SNR_values_db) max(SNR_values_db) 1/Nrealizations/10 1])
+else
+    disp(error)
+end
+
diff --git a/matlab/tests/pdcch_bler.m b/matlab/tests/pdcch_bler.m
new file mode 100644
index 000000000..c337a7b63
--- /dev/null
+++ b/matlab/tests/pdcch_bler.m
@@ -0,0 +1,171 @@
+%% PDCCH Blind Search and DCI Decoding + PCFICH encoding/decoding
+
+%% Cell-Wide Settings
+% A structure |enbConfig| is used to configure the eNodeB.
+clear
+
+Npackets = 1000;
+SNR_values = linspace(-0.5,3,6);
+
+txCFI = 3; 
+enbConfig.NDLRB = 15;                % No of Downlink RBs in total BW
+enbConfig.CyclicPrefix = 'Normal';  % CP length
+enbConfig.CFI = txCFI;                                                                                                                                                                                                                                                                                                                                                                                                                           ;                  % 4 PDCCH symbols as NDLRB <= 10
+enbConfig.Ng = 'Sixth';             % HICH groups
+enbConfig.CellRefP = 1;             % 1-antenna ports
+enbConfig.NCellID = 0;             % Physical layer cell identity
+enbConfig.NSubframe = 5;            % Subframe number 0
+enbConfig.DuplexMode = 'FDD';       % Frame structure
+enbConfig.PHICHDuration = 'Normal';
+C_RNTI = 65535;                         % 16-bit UE-specific mask
+
+%% Setup Fading channel model 
+cfg.Seed = 8;                  % Random channel seed
+cfg.NRxAnts = 1;               % 1 receive antenna
+cfg.DelayProfile = 'EPA';      % EVA delay spread
+cfg.DopplerFreq = 5;           % 120Hz Doppler frequency
+cfg.MIMOCorrelation = 'Low';   % Low (no) MIMO correlation
+cfg.InitTime = 0;              % Initialize at time zero
+cfg.NTerms = 16;               % Oscillators used in fading model
+cfg.ModelType = 'GMEDS';       % Rayleigh fading model type
+cfg.InitPhase = 'Random';      % Random initial phases
+cfg.NormalizePathGains = 'On'; % Normalize delay profile power 
+cfg.NormalizeTxAnts = 'On';    % Normalize for transmit antennas
+
+% Setup channel equalizer
+cec.PilotAverage = 'UserDefined';     % Type of pilot averaging
+cec.FreqWindow = 9;                   % Frequency window size
+cec.TimeWindow = 9;                   % Time window size
+cec.InterpType = 'linear';             % 2D interpolation type
+cec.InterpWindow = 'Centered';        % Interpolation window type
+cec.InterpWinSize = 1;                % Interpolation window size
+
+%% DCI Message Generation
+% Generate a DCI message to be mapped to the PDCCH.
+
+dciConfig.DCIFormat = 'Format1A';   % DCI message format
+dciConfig.Allocation.RIV = 26;      % Resource indication value
+
+% Create DCI message for given configuration
+[dciMessage, dciMessageBits] = lteDCI(enbConfig, dciConfig);
+
+%% DCI Channel Coding
+
+pdcchConfig.RNTI = C_RNTI;            % Radio network temporary identifier
+pdcchConfig.PDCCHFormat = 3;          % PDCCH format
+ueConfig.RNTI = C_RNTI;
+
+% DCI message bits coding to form coded DCI bits
+codedDciBits = lteDCIEncode(pdcchConfig, dciMessageBits);
+
+%% PDCCH Bits Generation
+
+pdcchDims = ltePDCCHInfo(enbConfig);
+
+% Initialize elements with -1 to indicate that all the bits are unused
+pdcchBits = -1*ones(pdcchDims.MTot, 1);
+
+% Perform search space for UE-specific control channel candidates.
+candidates = ltePDCCHSpace(enbConfig, pdcchConfig, {'bits', '1based'});
+
+Ncad=randi(length(candidates),1,1);
+
+% Map PDCCH payload on available UE-specific candidate. In this example the
+% first available candidate is used to map the coded DCI bits.
+pdcchBits ( candidates(Ncad, 1) : candidates(Ncad, 2) ) = codedDciBits;
+
+%% PDCCH Complex-Valued Modulated Symbol Generation
+
+pdcchSymbols = ltePDCCH(enbConfig, pdcchBits);
+pdcchIndices = ltePDCCHIndices(enbConfig,{'1based'});
+subframe_tx = lteDLResourceGrid(enbConfig);
+subframe_tx(pdcchIndices) = pdcchSymbols;
+
+%% PCFICH
+cfiCodeword = lteCFI(enbConfig);
+pcfichSymbols = ltePCFICH(enbConfig,cfiCodeword);
+pcfichIndices = ltePCFICHIndices(enbConfig,'1based');
+subframe_tx(pcfichIndices) = pcfichSymbols;
+
+%% Add references
+cellRsSym = lteCellRS(enbConfig);
+cellRsInd = lteCellRSIndices(enbConfig);
+subframe_tx(cellRsInd) = cellRsSym;
+
+[txWaveform, info] = lteOFDMModulate(enbConfig,subframe_tx);
+cfg.SamplingRate = info.SamplingRate;
+
+addpath('../../debug/lte/phy/lib/phch/test')
+
+decoded = zeros(size(SNR_values));
+decoded_cfi = zeros(size(SNR_values));
+decoded_liblte = zeros(size(SNR_values));
+decoded_cfi_liblte = zeros(size(SNR_values));
+
+parfor snr_idx=1:length(SNR_values)
+    SNRdB = SNR_values(snr_idx);
+    SNR = 10^(SNRdB/10);    % Linear SNR  
+    N0  = 1/(sqrt(2.0*enbConfig.CellRefP*double(info.Nfft))*SNR);
+    for i=1:Npackets
+
+        enbConfigRx=enbConfig; 
+
+        rxWaveform = sum(txWaveform,2);
+        
+        %% Fading
+        rxWaveform = lteFadingChannel(cfg,rxWaveform);
+        
+        %% Noise Addition
+        noise = N0*complex(randn(size(rxWaveform)), randn(size(rxWaveform)));  % Generate noise
+        rxWaveform = rxWaveform + noise; 
+       
+        %% Demodulate 
+        subframe_rx = lteOFDMDemodulate(enbConfigRx, rxWaveform);
+
+        % Perform channel estimation
+        [hest, nest] = lteDLChannelEstimate(enbConfigRx, cec, subframe_rx);
+
+        [pcfichSymbolsRx, pdcfichSymbolsHest] = lteExtractResources(pcfichIndices(:,1), subframe_rx, hest);
+
+        %% PCFICH decoding
+        [pcfichBits, pcfichSymbols] = ltePCFICHDecode(enbConfigRx,pcfichSymbolsRx, pdcfichSymbolsHest, nest);
+        rxCFI = lteCFIDecode(pcfichBits);
+        
+        decoded_cfi(snr_idx) = decoded_cfi(snr_idx) + (rxCFI == txCFI); 
+        
+        %% PDCCH Decoding
+        enbConfigRx.CFI = rxCFI; 
+        pdcchIndicesRx = ltePDCCHIndices(enbConfigRx,{'1based'});
+        [pdcchSymbolsRx, pdcchSymbolsHest] = lteExtractResources(pdcchIndicesRx(:,1), subframe_rx, hest);
+        [recPdcchBits] = ltePDCCHDecode(enbConfigRx, pdcchSymbolsRx, pdcchSymbolsHest, nest);
+
+        %% Blind Decoding using DCI Search        
+        [rxDCI, rxDCIBits] = ltePDCCHSearch(enbConfigRx, ueConfig, recPdcchBits);
+        decoded(snr_idx) = decoded(snr_idx) + (length(rxDCI)>0);
+
+        
+        %% Same with libLTE
+        [rxCFI, pcfichSymbols2, pcfichSymbolsRx2] = liblte_pcfich(enbConfigRx, rxWaveform);
+        decoded_cfi_liblte(snr_idx) = decoded_cfi_liblte(snr_idx) + (rxCFI == txCFI); 
+        enbConfigRx.CFI = rxCFI;
+        [found_liblte, llr, pdcchSymbols2] = liblte_pdcch(enbConfigRx, ueConfig.RNTI, rxWaveform);
+        decoded_liblte(snr_idx) = decoded_liblte(snr_idx)+found_liblte;
+    end
+    fprintf('SNR: %.1f\n',SNRdB)
+end
+
+if (Npackets>1)
+    semilogy(SNR_values,1-decoded/Npackets,'bo-',...
+             SNR_values,1-decoded_cfi/Npackets,'bx:',...
+             SNR_values,1-decoded_liblte/Npackets, 'ro-',...
+             SNR_values,1-decoded_cfi_liblte/Npackets,'rx:')
+    grid on
+    legend('Matlab all','Matlab cfi', 'libLTE all', 'libLTE cfi')
+    xlabel('SNR (dB)')
+    ylabel('BLER')
+    axis([min(SNR_values) max(SNR_values) 1/Npackets/10 1])
+else
+    disp(decoded)
+    disp(decoded_liblte)
+end
+
diff --git a/matlab/tests/pdsch_resalloc.m b/matlab/tests/pdsch_resalloc.m
new file mode 100644
index 000000000..61c9d9740
--- /dev/null
+++ b/matlab/tests/pdsch_resalloc.m
@@ -0,0 +1,76 @@
+filename='../../debug/dist_ra.dat';
+
+enb.NDLRB = 50;
+enb.CyclicPrefix = 'Normal';
+enb.PHICHDuration = 'Normal';
+enb.CFI = 2;
+enb.Ng = 'Sixth';
+enb.CellRefP = 1;
+enb.NCellID = 196;
+enb.NSubframe = 5;
+enb.NTotalSubframes=1;
+enb.DuplexMode = 'FDD';
+
+dci.NDLRB = enb.NDLRB;
+dci.DCIFormat = 'Format1C';
+dci.AllocationType=1;
+%dci.Allocation.Bitmap='01111000011110000';
+%dci.Allocation.Subset=3;
+dci.Allocation.RIV = 33;
+dci.Allocation.Gap = 0;
+dci.ModCoding=6;
+dci.RV=0;
+dci.DuplexMode = enb.DuplexMode;
+dci.NTxAnts = enb.CellRefP;
+pdcch.RNTI = 65535;
+pdcch.PDCCHFormat = 3;
+
+pdsch.Modulation='QPSK';
+pdsch.RNTI=pdcch.RNTI;
+if (enb.CellRefP == 1) 
+    pdsch.TxScheme='Port0';
+else
+    pdsch.TxScheme='TxDiversity';
+end
+pdsch.NLayers=enb.CellRefP;
+pdsch.trblklen=176;
+pdsch.RV=dci.RV;
+
+% Begin frame generation
+subframe = lteDLResourceGrid(enb);
+
+%%% Create Reference Signals
+rsAnt = lteCellRS(enb);
+indAnt = lteCellRSIndices(enb);
+subframe(indAnt) = rsAnt;
+
+%%% Create PDCCH
+[dciMessage,dciMessageBits] = lteDCI(enb,dci);
+codedDciBits = lteDCIEncode(pdcch,dciMessageBits);
+pdcchInfo = ltePDCCHInfo(enb);
+pdcchBits = -1*ones(1,pdcchInfo.MTot);
+candidates = ltePDCCHSpace(enb,pdcch,{'bits','1based'});
+pdcchBits (candidates(1,1):candidates(1,2)) = codedDciBits;
+pdcchSymbols = ltePDCCH(enb, pdcchBits);
+pdcchIndices = ltePDCCHIndices(enb,{'1based'});
+subframe(pdcchIndices) = pdcchSymbols;
+
+% Create PDSCH
+pdsch.prbset = lteDCIResourceAllocation(enb,dci);
+
+[pdschIndices,pdschInfo] = ltePDSCHIndices(enb,pdsch,pdsch.prbset);
+
+dlschTransportBlk=randi([0 1],pdsch.trblklen,1);
+pdschcodeword = lteDLSCH(enb,pdsch,pdschInfo.G,dlschTransportBlk);
+%crced = lteCRCEncode(dlschTransportBlk, '24A');
+%encoded = lteTurboEncode(crced);
+%pdschcodeword2 = lteRateMatchTurbo(encoded,pdschInfo.G,pdsch.RV); 
+pdschSymbols = ltePDSCH(enb,pdsch,pdschcodeword);
+
+subframe(pdschIndices) = pdschSymbols;
+
+txwaveform = lteOFDMModulate(enb,subframe);
+
+write_complex(filename,sum(txwaveform,2));
+fprintf('Written signal to %s\n',filename);
+
diff --git a/matlab/tests/sss_test.m b/matlab/tests/sss_test.m
new file mode 100644
index 000000000..21c9a0bac
--- /dev/null
+++ b/matlab/tests/sss_test.m
@@ -0,0 +1,144 @@
+
+SNR_values = linspace(-6,4,10);
+Npackets = 200; 
+CFO=4/15;
+m0=7;
+m1=10;
+%m0=26;
+%m1=21;
+
+recordedWaveform = x;
+if (~isempty(recordedWaveform))
+    Npackets = floor(length(recordedWaveform)/19200)-1;
+    SNR_values = 0;
+end
+
+error = zeros(6,length(SNR_values));
+
+enb = struct('NCellID',2,'NSubframe',0,'NDLRB',6,'CellRefP',1,'CyclicPrefix','Normal','DuplexMode','FDD');
+sss=lteSSS(enb);
+
+cfg.Seed = 2;                  % Random channel seed
+cfg.NRxAnts = 1;               % 1 receive antenna
+cfg.DelayProfile = 'ETU';      % EVA delay spread
+cfg.DopplerFreq = 144;         % 120Hz Doppler frequency
+cfg.MIMOCorrelation = 'Low';   % Low (no) MIMO correlation
+cfg.NTerms = 16;               % Oscillators used in fading model
+cfg.ModelType = 'GMEDS';       % Rayleigh fading model type
+cfg.InitPhase = 'Random';      % Random initial phases
+cfg.NormalizePathGains = 'On'; % Normalize delay profile power 
+cfg.NormalizeTxAnts = 'On';    % Normalize for transmit antennas       % Initialize at time zero
+
+[s, c0, c1] = get_sc(mod(enb.NCellID,3));
+
+subframe = lteDLResourceGrid(enb);
+sssSym = lteSSS(enb);
+sssInd = lteSSSIndices(enb);
+subframe(sssInd) = sssSym;
+N_id_1 = floor(enb.NCellID/3);
+
+[txWaveform,info] = lteOFDMModulate(enb,subframe); 
+cfg.SamplingRate = info.SamplingRate;
+fftSize = info.Nfft;
+
+
+addpath('../../debug/lte/phy/lib/sync/test')
+
+for snr_idx=1:length(SNR_values)
+    SNRdB = SNR_values(snr_idx);
+    for i=1:Npackets
+        %% Noise Addition
+        SNR = 10^(SNRdB/10);    % Linear SNR  
+
+        if (isempty(recordedWaveform))
+            cfg.InitTime = i*(10^-3);              
+            [rxWaveform, info]= lteFadingChannel(cfg,txWaveform);
+            rxWaveform = txWaveform; 
+
+            % Add CFO
+            freq = CFO/double(fftSize);
+            rxWaveform = rxWaveform.*exp(1i*2*pi*freq*(1:length(txWaveform))');
+
+            N0 = 1/(sqrt(2.0*enb.CellRefP*double(fftSize))*SNR);
+            noise = N0*complex(randn(size(rxWaveform)), randn(size(rxWaveform)));  % Generate noise
+
+            rxWaveform = rxWaveform + noise;
+        else
+            rxWaveform = recordedWaveform(i*19200+1:(i+1)*19200); 
+        end
+        
+        offset = lteDLFrameOffset(enb,rxWaveform);
+        offsetVec(i)=offset;
+        rxWaveform = [rxWaveform(1+offset:end,:); zeros(offset,1)];
+
+        subframe_rx = lteOFDMDemodulate(enb,rxWaveform,1);
+        
+        sss_rx = subframe_rx(lteSSSIndices(enb));
+        sss0=sss_rx(1:2:end);
+        sss1=sss_rx(2:2:end);
+
+        beta0=sss0.*c0';
+        beta1=sss1.*c1';
+
+        corr0=zeros(31,1);
+        for m=1:31
+            corr0(m)=sum(beta0.*s(m,:)');    
+        end
+        corr0=abs(corr0).^2;
+        [m, idx]=max(corr0);
+        
+        error(1,snr_idx) = error(1,snr_idx) + ((idx ~= m0 && idx ~= m1));
+        
+        M=2;
+        Nm=10;
+
+        corr2=zeros(31,1);
+        for m=1:31
+            for j=0:M
+                idx=1+j*Nm:(j+1)*Nm;
+                corr2(m)=corr2(m)+abs(sum(beta0(idx).*s(m,idx)')).^2;
+            end
+        end
+        [m, idx]=max(corr2);
+        
+        error(2,snr_idx) = error(2,snr_idx) + ((idx ~= m0 && idx ~= m1));
+
+        corr3=zeros(31,1);
+        for m=1:31
+            corr3(m)=abs(sum(beta0(2:end).*conj(beta0(1:end-1)).*transpose(s(m,2:end).*conj(s(m,1:end-1))))).^2;            
+        end
+        [m, idx]=max(corr3);
+        
+        error(3,snr_idx) = error(3,snr_idx) + ((idx ~= m0 && idx ~= m1));
+        
+        % libLTE results
+        [n,sf_idx,lt_corr0]=liblte_sss(enb,rxWaveform,'full');
+        [m, idx]=max(lt_corr0);
+        error(4,snr_idx) = error(4,snr_idx) + ((idx ~= m0 && idx ~= m1));
+
+        [n,sf_idx,lt_corr2]=liblte_sss(enb,rxWaveform,'partial');
+        [m, idx]=max(lt_corr2);
+        error(5,snr_idx) = error(5,snr_idx) + ((idx ~= m0 && idx ~= m1));
+        
+        [n,sf_idx,lt_corr3]=liblte_sss(enb,rxWaveform,'diff');
+        [m, idx]=max(lt_corr3);
+        error(6,snr_idx) = error(6,snr_idx) + ((idx ~= m0 && idx ~= m1));
+        
+    end
+end
+
+if (length(SNR_values) > 1) 
+    plot(SNR_values,1-error/Npackets)
+    legend('Full','Partial','Differential','Full-lt','Partial-lt','Differential-lt')
+    grid on
+else
+    e=error/Npackets;
+    fprintf('Full (mt/lt): \t%f/%f\n',e(1),e(4));
+    fprintf('Partial (mt/lt):%f/%f\n',e(2),e(5));
+    fprintf('Diff (mt/lt): \t%f/%f\n',e(3),e(6));
+end
+
+
+
+
+
diff --git a/matlab/tests/viterbi_bler.m b/matlab/tests/viterbi_bler.m
new file mode 100644
index 000000000..06b49c82f
--- /dev/null
+++ b/matlab/tests/viterbi_bler.m
@@ -0,0 +1,49 @@
+
+clear
+blen=40;
+SNR_values_db=linspace(-6,4,8);
+Nrealizations=5000;
+
+addpath('../../debug/lte/phy/lib/fec/test')
+
+errors1=zeros(1,length(SNR_values_db));
+errors2=zeros(1,length(SNR_values_db));
+for snr_idx=1:length(SNR_values_db)
+    SNRdB = SNR_values_db(snr_idx);             % Desired SNR in dB
+    SNR = 10^(SNRdB/20);    % Linear SNR  
+
+    for i=1:Nrealizations
+        Data = randi(2,blen,1)==1;
+        codedData = lteConvolutionalEncode(Data);
+
+        codedsymbols = 2*double(codedData)-1;
+        
+        %% Additive Noise
+        N0 = 1/SNR;
+
+        % Create additive white Gaussian noise
+        noise = N0*randn(size(codedsymbols));   
+        
+        noisysymbols = noise + codedsymbols; 
+        
+        decodedData = lteConvolutionalDecode(noisysymbols);
+        interleavedSymbols = reshape(reshape(noisysymbols,[],3)',1,[]);
+        [decodedData2, quant] = liblte_viterbi(interleavedSymbols);
+        
+        errors1(snr_idx) = errors1(snr_idx) + any(decodedData ~= Data);
+        errors2(snr_idx) = errors2(snr_idx) + any(decodedData2 ~= Data);
+    end
+end
+
+if (length(SNR_values_db) > 1)
+    semilogy(SNR_values_db, errors1/Nrealizations, ...
+    SNR_values_db, errors2/Nrealizations)
+    grid on
+    xlabel('SNR (dB)')
+    ylabel('BLER')
+    legend('Matlab','libLTE');
+else
+    disp(errors1);
+    disp(errors2);
+    disp(errors3);
+end
\ No newline at end of file