Gravitational waves from oscillating neutron stars in axial symmetry are studied performing numerical simulations in full general relativity. Neutron stars are modeled by a polytropic equation of state for simplicity. A gauge-invariant wave extraction method as well as a quadrupole formula are adopted for computation of gravitational waves. It is found that the gauge-invariant variables systematically contain numerical errors generated near the outer boundaries in the present axisymmetric computation. We clarify their origin, and illustrate that it is possible to eliminate the dominant part of the systematic errors. The best corrected waveforms for oscillating and rotating stars currently contain errors of magnitude similar to10(-3) in the local wave zone. Comparing the waveforms obtained by the gauge-invariant technique with those by the quadrupole formula, it is shown that the quadrupole formula yields approximate gravitational waveforms in addition to a systematic underestimation of the amplitude of O(M/R) where M and R denote the mass and the radius of neutron stars. However, the wave phase and modulation of the amplitude can be computed accurately. This indicates that the quadrupole formula is a useful tool for studying gravitational waves from rotating stellar core collapse to a neutron star in fully general relativistic simulations. The properties of the gravitational waveforms from the oscillating and rigidly rotating neutron stars are also addressed paying attention to the oscillation associated with fundamental modes.