{"cells":[{"metadata":{},"cell_type":"markdown","source":"This notebook compares the speed of three methods for resampling audio data: librosa, the existing torchaudio resample function and a new, efficient pytorch resampler (code below) based on https://github.com/danpovey/filtering/blob/master/lilfilter/resampler.py\n\n\nComparison 1: 'kaiser_best' settings in librosa vs 'kaiser_best' setting in the efficient pytorch resampler (should be the same set-up)\n* librosa: 51 s \n* efficient pytorch resampler: 9 s\n\nComparison 2: default setting in torchaudio vs window='hann', num_zeros=6 in the efficient pytorch resampler (should be the same set-up)\n* torchaudio: 10 s\n* efficient pytorch resampler: 1 s\n"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nimport matplotlib.pyplot as plt\nimport librosa\nimport time\nimport torchaudio\nimport torch\n\nimport warnings\nwarnings.filterwarnings(\"ignore\", message=\"PySoundFile failed. Trying audioread instead.\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import torch\nimport numpy as np\nfrom scipy import special\n\n\nclass Resampler(torch.nn.Module):\n    \"\"\"\n    Efficiently resample audio signals\n    This module is much faster than resampling with librosa because it exploits pytorch's efficient conv1d operations\n    This module is also faster than the existing pytorch resample function in\n    https://github.com/pytorch/audio/blob/b6a61c3f7d0267c77f8626167cc1eda0335f2753/torchaudio/compliance/kaldi.py#L892\n    \n    Based on \n    https://github.com/danpovey/filtering/blob/master/lilfilter/resampler.py\n    with improvements to include additional filter types and input parameters that align with the librosa api\n    \"\"\"\n\n    def __init__(self,\n                 input_sr, output_sr, dtype,\n                 num_zeros = 64, cutoff_ratio = 0.95, filter='kaiser', beta=14.0):\n        super().__init__()  # init the base class\n        \"\"\"\n        This creates an object that can apply a symmetric FIR filter\n        based on torch.nn.functional.conv1d.\n\n        Args:\n          input_sr:  The input sampling rate, AS AN INTEGER..\n              does not have to be the real sampling rate but should\n              have the correct ratio with output_sr.\n          output_sr:  The output sampling rate, AS AN INTEGER.\n              It is the ratio with the input sampling rate that is\n              important here.\n          dtype:  The torch dtype to use for computations (would be preferrable to \n               set things up so passing the dtype isn't necessary)\n          num_zeros: The number of zeros per side in the (sinc*hanning-window)\n              filter function.  More is more accurate, but 64 is already\n              quite a lot. The kernel size is 2*num_zeros + 1.\n          cutoff_ratio: The filter rolloff point as a fraction of the\n             Nyquist frequency.\n          filter: one of ['kaiser', 'kaiser_best', 'kaiser_fast', 'hann']\n          beta: parameter for 'kaiser' filter\n\n        You can think of this algorithm as dividing up the signals\n        (input,output) into blocks where there are `input_sr` input\n        samples and `output_sr` output samples.  Then we treat it\n        using convolutional code, imagining there are `input_sr`\n        input channels and `output_sr` output channels per time step.\n\n        \"\"\"\n        assert isinstance(input_sr, int) and isinstance(output_sr, int)\n        if input_sr == output_sr:\n            self.resample_type = 'trivial'\n            return\n        \n        def gcd(a, b):\n            \"\"\" Return the greatest common divisor of a and b\"\"\"\n            assert isinstance(a, int) and isinstance(b, int)\n            if b == 0:\n                return a\n            else:\n                return gcd(b, a % b)  \n\n        d = gcd(input_sr, output_sr)\n        input_sr, output_sr = input_sr // d, output_sr // d\n\n        assert dtype in [torch.float32, torch.float64]\n        assert num_zeros > 3  # a reasonable bare minimum\n        np_dtype = np.float32 if dtype == torch.float32 else np.float64\n\n        assert filter in ['hann', 'kaiser', 'kaiser_best', 'kaiser_fast']\n\n        if filter == 'kaiser_best':\n            num_zeros = 64\n            beta = 14.769656459379492\n            cutoff_ratio = 0.9475937167399596\n            filter = 'kaiser'\n        elif filter == 'kaiser_fast':\n            num_zeros = 16\n            beta = 8.555504641634386\n            cutoff_ratio = 0.85\n            filter = 'kaiser'\n\n        # Define one 'block' of samples `input_sr` input samples\n        # and `output_sr` output samples.  We can divide up\n        # the samples into these blocks and have the blocks be\n        #in correspondence.\n\n        # The sinc function will have, on average, `zeros_per_block`\n        # zeros per block.\n        zeros_per_block = min(input_sr, output_sr) * cutoff_ratio\n\n        # The convolutional kernel size will be n = (blocks_per_side*2 + 1),\n        # i.e. we add that many blocks on each side of the central block.  The\n        # window radius (defined as distance from center to edge)\n        # is `blocks_per_side` blocks.  This ensures that each sample in the\n        # central block can \"see\" all the samples in its window.\n        #\n        # Assuming the following division is not exact, adding 1\n        # will have the same effect as rounding up.\n        #blocks_per_side = 1 + int(num_zeros / zeros_per_block)\n        blocks_per_side = int(np.ceil(num_zeros / zeros_per_block))\n\n        kernel_width = 2*blocks_per_side + 1\n\n        # We want the weights as used by torch's conv1d code; format is\n        #  (out_channels, in_channels, kernel_width)\n        # https://pytorch.org/docs/stable/nn.functional.html\n        weights = torch.tensor((output_sr, input_sr, kernel_width), dtype=dtype)\n\n        # Computations involving time will be in units of 1 block.  Actually this\n        # is the same as the `canonical` time axis since each block has input_sr\n        # input samples, so it would be one of whatever time unit we are using\n        window_radius_in_blocks = blocks_per_side\n\n\n        # The `times` below will end up being the args to the sinc function.\n        # For the shapes of the things below, look at the args to `view`.  The terms\n        # below will get expanded to shape (output_sr, input_sr, kernel_width) through\n        # broadcasting\n        # We want it so that, assuming input_sr == output_sr, along the diagonal of\n        # the central block we have t == 0.\n        # The signs of the output_sr and input_sr terms need to be opposite.  The\n        # sign that the kernel_width term needs to be will depend on whether it's\n        # convolution or correlation, and the logic is tricky.. I will just find\n        # which sign works.\n\n\n        times = (\n            np.arange(output_sr, dtype=np_dtype).reshape((output_sr, 1, 1)) / output_sr -\n            np.arange(input_sr, dtype=np_dtype).reshape((1, input_sr, 1)) / input_sr -\n            (np.arange(kernel_width, dtype=np_dtype).reshape((1, 1, kernel_width)) - blocks_per_side))\n\n\n        def hann_window(a):\n            \"\"\"\n            hann_window returns the Hann window on [-1,1], which is zero\n            if a < -1 or a > 1, and otherwise 0.5 + 0.5 cos(a*pi).\n            This is applied elementwise to a, which should be a NumPy array.\n\n            The heaviside function returns (a > 0 ? 1 : 0).\n            \"\"\"\n            return np.heaviside(1 - np.abs(a), 0.0) * (0.5 + 0.5 * np.cos(a * np.pi))\n\n        def kaiser_window(a, beta):\n            w = special.i0(beta * np.sqrt(np.clip(1 - ((a - 0.0) / 1.0) ** 2.0, 0.0, 1.0))) / special.i0(beta)\n            return np.heaviside(1 - np.abs(a), 0.0) * w\n\n\n        # The weights below are a sinc function times a Hann-window function.\n        #\n        # Multiplication by zeros_per_block normalizes the sinc function \n        # (to compensate for scaling on the x-axis), so that the integral is 1.\n        #\n        # Division by input_sr normalizes the input function. Think of the input \n        # as a stream of dirac deltas passing through a low pass filter: \n        # in order to have the same magnitude as the original input function, \n        # we need to divide by the number of those deltas per unit time.\n        if filter == 'hann':\n            weights = (np.sinc(times * zeros_per_block)\n                       * hann_window(times / window_radius_in_blocks)\n                       * zeros_per_block / input_sr)\n        else:\n            weights = (np.sinc(times * zeros_per_block)\n                       * kaiser_window(times / window_radius_in_blocks, beta)\n                       * zeros_per_block / input_sr)\n\n        self.input_sr = input_sr\n        self.output_sr = output_sr\n\n        # weights has dim (output_sr, input_sr, kernel_width).  \n        # If output_sr == 1, we can fold the input_sr into the\n        # kernel_width (i.e. have just 1 input channel); this will make the\n        # convolution faster and avoid unnecessary reshaping.\n\n        assert weights.shape == (output_sr, input_sr, kernel_width)\n        if output_sr == 1:\n            self.resample_type = 'integer_downsample'\n            self.padding = input_sr * blocks_per_side\n            weights = torch.tensor(weights, dtype=dtype, requires_grad=False)\n            self.weights = weights.transpose(1, 2).contiguous().view(1, 1, input_sr * kernel_width)\n\n        elif input_sr == 1:\n            # In this case we'll be doing conv_transpose, so we want the same weights that\n            # we would have if we were *downsampling* by this factor-- i.e. as if input_sr,\n            # output_sr had been swapped.\n            self.resample_type = 'integer_upsample'\n            self.padding = output_sr * blocks_per_side\n            weights = torch.tensor(weights, dtype=dtype, requires_grad=False)\n            self.weights = weights.flip(2).transpose(0, 2).contiguous().view(1, 1, output_sr * kernel_width)\n        else:\n            self.resample_type = 'general'\n            self.reshaped = False\n            self.padding = blocks_per_side\n            self.weights = torch.tensor(weights, dtype=dtype, requires_grad=False)\n\n        self.weights = torch.nn.Parameter(self.weights, requires_grad=False)      \n\n    @torch.no_grad()\n    def forward(self, data):\n        \"\"\"\n        Resample the data\n\n        Args:\n         input: a torch.Tensor with the same dtype as was passed to the\n           constructor.\n         There must be 2 axes, interpreted as (minibatch_size, sequence_length)...\n         the minibatch_size may in practice be the number of channels.\n\n        Return:  Returns a torch.Tensor with the same dtype as the input, and\n         dimension (minibatch_size, (sequence_length//input_sr)*output_sr),\n         where input_sr and output_sr are the corresponding constructor args,\n         modified to remove any common factors.\n        \"\"\"\n        if self.resample_type == 'trivial':\n            return data\n        elif self.resample_type == 'integer_downsample':\n            (minibatch_size, seq_len) = data.shape\n            # will be shape (minibatch_size, in_channels, seq_len) with in_channels == 1\n            data = data.unsqueeze(1)\n            data = torch.nn.functional.conv1d(data,\n                                             self.weights,\n                                             stride=self.input_sr,\n                                             padding=self.padding)\n            # shape will be (minibatch_size, out_channels = 1, seq_len);\n            # return as (minibatch_size, seq_len)\n            return data.squeeze(1)\n\n        elif self.resample_type == 'integer_upsample':\n            data = data.unsqueeze(1)\n            data = torch.nn.functional.conv_transpose1d(data,\n                                                      self.weights,\n                                                      stride=self.output_sr,\n                                                      padding=self.padding)\n\n            return data.squeeze(1)\n        else:\n            assert self.resample_type == 'general'\n            (minibatch_size, seq_len) = data.shape\n            num_blocks = seq_len // self.input_sr\n            if num_blocks == 0:\n                # TODO: pad with zeros.\n                raise RuntimeError(\"Signal is too short to resample\")\n            #data = data[:, 0:(num_blocks*self.input_sr)]  # Truncate input\n            data = data[:, 0:(num_blocks*self.input_sr)].view(minibatch_size, num_blocks, self.input_sr)\n\n\n            # Torch's conv1d expects input data with shape (minibatch, in_channels, time_steps), so transpose\n            data = data.transpose(1, 2)\n\n\n            data = torch.nn.functional.conv1d(data, self.weights,\n                                          padding=self.padding)\n\n            assert data.shape == (minibatch_size, self.output_sr, num_blocks)\n            return data.transpose(1, 2).contiguous().view(minibatch_size, num_blocks * self.output_sr)\n\n\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"base_path = \"../input/birdsong-recognition/\"\ntrain_df = pd.read_csv(os.path.join(base_path, \"train.csv\"))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#convert the sample rates to int from string\ntrain_df['sampling_rate'] = train_df['sampling_rate'].apply(lambda x: int(x.split(' ')[0]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"'''\nsample rate census for the first 20 files\n'''\ntrain_df.loc[0:19, 'sampling_rate'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"'''\nload the first 20 audio files\n'''\n\naudio = {}\nsamplerates = {}\nfor i in range(20):\n    fname = base_path + \"train_audio/\" + train_df[\"ebird_code\"].iloc[i] +\"/\"+ train_df[\"filename\"].iloc[i]\n    y, sr = librosa.load(fname, sr=None)\n    audio[i] = y\n    samplerates[i] = sr","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"'''\nresample with librosa, 'kaiser_best' settings\n'''\n\nt0 = time.time()\nfor i in range(20):\n    y_librosa = librosa.resample(audio[i], samplerates[i], 22050, res_type='kaiser_best')\nt1 = time.time()\n\nprint(f'execution time librosa resample (kaiser_best): {t1-t0}')    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"'''\nresample with efficient pytorch resampler\n'kaiser_best' reproduces the 'kaiser_best' settings in librosa\n'''\n\nresampler = {}\nresampler[44100] = Resampler(input_sr=44100, output_sr=22050, dtype=torch.float32, filter='kaiser_best')\nresampler[48000] = Resampler(input_sr=48000, output_sr=22050, dtype=torch.float32, filter='kaiser_best')\n\nt0 = time.time()\nfor i in range(20):\n    if len(audio[i].shape) == 1:\n        y = torch.tensor(audio[i]).unsqueeze(0)\n    else:\n        y = torch.tensor(audio[i])\n   \n    y_best = resampler[samplerates[i]].forward(y)\nt1 = time.time()\n\nprint(f'execution time efficient torch resample (kaiser_best): {t1-t0}') ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"'''\nresample with torchaudio's resample\n'''\n\nresampler = {}\nresampler[44100] = torchaudio.transforms.Resample(orig_freq=44100, new_freq=22050)\nresampler[48000] = torchaudio.transforms.Resample(orig_freq=48000, new_freq=22050)\n\nt0 = time.time()\nfor i in range(20):\n    y_torchaudio = resampler[samplerates[i]].forward(torch.tensor(audio[i]))\nt1 = time.time()\n\nprint(f'execution time torchaudio resample: {t1-t0}') ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"'''\nresample with efficient pytorch resampler\nfilter='hann' and num_zeros=6 reproduce the settings in torchaudio's existing resample function from above\n'''\n\nresampler = {}\nresampler[44100] = Resampler(input_sr=44100, output_sr=22050, dtype=torch.float32, filter='hann', num_zeros=6)\nresampler[48000] = Resampler(input_sr=48000, output_sr=22050, dtype=torch.float32, filter='hann', num_zeros=6)\n\nt0 = time.time()\nfor i in range(20):\n    if len(audio[i].shape) == 1:\n        y = torch.tensor(audio[i]).unsqueeze(0)\n    else:\n        y = torch.tensor(audio[i])\n  \n    y_best_hann = resampler[samplerates[i]].forward(y)\nt1 = time.time()\n\nprint(f'execution time efficient torch resample (hann, num_zeros=6): {t1-t0}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(y_librosa.shape)\nprint(y_best.size())\nprint(y_torchaudio.size())\nprint(y_best_hann.size())\n\n\ntorch.tensor(y_librosa)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(np.max(y_librosa))\nprint(np.max(y_best.numpy()))\nprint(np.max(np.abs(y_librosa - y_best.squeeze(0).numpy())))\nprint(np.mean(np.abs(y_librosa - y_best.squeeze(0).numpy())))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(np.arange(len(y_librosa)), y_librosa, 'k.')\nplt.plot(np.arange(len(y_librosa)), y_best.squeeze(0).numpy(), 'r.')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(y_best.squeeze(0).numpy(), y_best.squeeze(0).numpy()-y_librosa, 'k.')\nplt.xlabel('y value efficient resample')\nplt.ylabel('y value (efficient resample - librosa resample)')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(np.max(y_torchaudio.squeeze(0).numpy()))\nprint(np.max(y_best_hann.numpy()))\nprint(np.max(np.abs(y_torchaudio.squeeze(0).numpy() - y_best_hann.squeeze(0).numpy())))\nprint(np.mean(np.abs(y_torchaudio.squeeze(0).numpy() - y_best_hann.squeeze(0).numpy())))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(y_best_hann.squeeze(0).numpy(), y_best_hann.squeeze(0).numpy()-y_torchaudio.squeeze(0).numpy(), 'k.')\nplt.xlabel('y value efficient resample')\nplt.ylabel('y value (efficient resample - torchaudio resample)')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"torch.sqrt(torch.mean(torch.pow(y_torchaudio, 2)))","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}