<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def demodulate(self, f, stride=1, exp=False, deg=True):
 """Compute the average magnitude and phase of this `TimeSeries` once per stride at a given frequency. Parameters f : `float` frequency (Hz) at which to demodulate the signal stride : `float`, optional stride (seconds) between calculations, defaults to 1 second exp : `bool`, optional return the magnitude and phase trends as one `TimeSeries` object representing a complex exponential, default: False deg : `bool`, optional if `exp=False`, calculates the phase in degrees Returns ------- mag, phase : `TimeSeries` if `exp=False`, returns a pair of `TimeSeries` objects representing magnitude and phase trends with `dt=stride` out : `TimeSeries` if `exp=True`, returns a single `TimeSeries` with magnitude and phase trends represented as `mag * exp(1j*phase)` with `dt=stride` Examples -------- Demodulation is useful when trying to examine steady sinusoidal signals we know to be contained within data. For instance, we can download some data from LOSC to look at trends of the amplitude and phase of LIGO Livingston's calibration line at 331.3 Hz: We can demodulate the `TimeSeries` at 331.3 Hz with a stride of one minute: We can then plot these trends to visualize fluctuations in the amplitude of the calibration line: """

        stridesamp = int(stride * self.sample_rate.value)
        nsteps = int(self.size // stridesamp)
        # stride through the TimeSeries and mix with a local oscillator,
        # taking the average over each stride
        out = type(self)(numpy.zeros(nsteps, dtype=complex))
        out.__array_finalize__(self)
        out.sample_rate = 1 / float(stride)
        w = 2 * numpy.pi * f * self.dt.decompose().value
        for step in range(nsteps):
            istart = int(stridesamp * step)
            iend = istart + stridesamp
            idx = numpy.arange(istart, iend)
            mixed = 2 * numpy.exp(-1j * w * idx) * self.value[idx]
            out.value[step] = mixed.mean()
        if exp:
            return out
        mag = out.abs()
        phase = type(mag)(numpy.angle(out, deg=deg))
        phase.__array_finalize__(out)
        phase.override_unit('deg' if deg else 'rad')
        return (mag, phase)