The design of informatively rich input signals is essential for accurate system identification, yet classical Fisher-information-based methods are inherently local and often inadequate in the presence of significant model uncertainty and nonlinearity. This paper develops a Bayesian approach that uses the mutual information (MI) between observations and parameters as the utility function. To address the computational intractability of the MI, we maximize a tractable MI lower bound. The method is then applied to the design of an input signals for the identification of quasi-linear stochastic dynamical systems. Evaluating the MI lower bound requires inversion of large covariance matrices whose dimensions scale with the number of data points $N$. To overcome this problem, an algorithm that reduces the dimension of the matrices to be inverted by a factor of $N$ is developed, making the approach feasible for long experiments. The proposed Bayesian method is compared with the average D-optimal design method, a semi-Bayesian approach, and its advantages are demonstrated. The effectiveness of the proposed method is further illustrated through four examples, including atomic sensor models, where the input signals that generates large MI are especially important for reducing the estimation error.