We develop a non-parametric method to reconstruct the primordial power spectrum PR(k) from Large Scale Structure (LSS) data through Bayesian inference and nested sampling. The performance of the method is studied by applying it to simulations of the clustering of two different object catalogues, low-z (ELGs) and high-z (QSOs), with two different photometric errors. These object clusterings are derived from different templates of the primordial power spectrum motivated by models of inflation: the Standard Model (SM) power law characterized by the two parameters As and ns; a local feature template; and a global oscillatory template. Our reconstruction method involves sampling N knots in the log {k,PR(k)} plane. We use two statistical tests to examine the reconstructions for signs of primordial features: a global test comparing the evidences and a novel local test quantifying the power of the hypothesis test between the power law model and the marginalized probability over N-knots model. Scenarios with different redshift bins, photometric errors, feature amplitudes and detection levels are discussed. The method shows good performance in all scenarios considered. In particular, the tests show no feature detection for the SM. The method is able to detect deviations from the power law power spectrum for all considered features with amplitude above a given threshold, and the best sensitivity is reached when combining several redshift bins of low-z or high-z objects. In addition, we include a first application to real data from the Sloan Digital Sky Survey Luminous Red Galaxy Data Release 4 (SDSS LRG 04), finding no preference for deviations from the primordial power law. The method is flexible and model independent, and progress on varying the cosmological parameters is being made to improve the completeness of the approach. Applications and forecasts for upcoming LSS catalogues, as Euclid, BOSS or J-PAS, are expected.