What you seem to want is a multiplicative modification of the spectrum of the input signal in the following form:
where $X(\omega)$ is the spectrum of the input signal, $Y(\omega)$ is the spectrum of the output signal, and $H(\omega)$ is the function that modifies the input spectrum.
Equation $(1)$ is exactly what a (linear and time-invariant) filter does. Examples of such filters in the discrete domain are FIR and IIR filters. In general, the function $H(\omega)$ - called the filter's frequency response - is complex-valued. That means that not only the magnitude of the input spectrum is changed but also its phase. If I understand your question correctly, you just care about the magnitude.
For clarification, the frequency response's influence on a sinusoidal input signal $x[n]=\sin(\omega_0n)$ is described by the following equation:
where $\phi(\omega)$ is the phase of the complex-valued frequency response $H(\omega)$:
Eq. $(2)$ shows that the magnitude of a sinusoid at frequency $\omega_0$ is modified by the magnitude of the frequency response at that frequency. Furthermore, its phase is modified by the phase of the frequency response at that frequency.
The catch is that you can't just prescribe some $H(\omega)$ and expect a realizable filter to exactly implement the given response. You will always have to live with an approximation of the desired response. However, if computational cost and memory requirements are of no concern you can get very close to the prescribed response.
I would recommend to use an FIR filter to approximate the given desired response. Reasons why one might want to avoid FIR filters are their large delay (at least for linear phase filters) and their large computational complexity, both when compared to equivalent IIR filters. However, I think that for your application these drawbacks are not very relevant, and the advantages of FIR filters outweigh the disadvantages. The advantages are ease of design and inherent stability.
The frequency response of a causal FIR filter is
where $h[n]$ are the filter coefficients ("taps"), and $N$ is the filter length (= number of taps). Obviously, the larger $N$, i.e., the more taps, the better the approximation of a given response will be.
Given the filter coefficients $h[n]$, the output $y[n]$ is computed by the discrete convolution of the filter coefficients with the input $x[n]$:
Designing a filter means finding the "best" filter coefficients $h[n]$ for a given specification. There are countless filter design methods and it's not always straightforward to find the one that best fits a given problem. In your case I would probably start experimenting with the frequency sampling method. In this case, the desired response is sampled on a dense equidistant frequency grid. This discretized response is then transformed to the time domain by applying an inverse discrete Fourier transform (IDFT), which is usually implemented by the (inverse) Fast Fourier Transform (IFFT). Finally, the resulting impulse response is usually windowed to smooth out the resulting response and to reduce the impulse response to the desired filter length $N$.
For illustration purposes I used the frequency sampling method to design a filter approximating a rather silly fantasy specification. The desired magnitude response has a staircase shape and was sampled at $1000$ equidistant points. The desired filter length was chosen to be $N=201$. The figure below shows the result. Clearly, an actual filter cannot perfectly follow the discontinuities of the desired response. However, increasing the filter length will make the transitions narrower.
is a link to a simple Octave/Matlab implementation of the frequency sampling method used in the design example shown above. The filter was designed with the following commands:
mag = [ones(1,200),.75*ones(1,200),.5*ones(1,200),.25*ones(1,200),zeros(1,200)];
h = fsamp( mag, 201 );