Bottom Line
If the OP intended to "down-sample", that is not at all what is being done. The OP is reducing the number of samples in the frequency domain. This is not time domain down-sampling or decimation. If this was the intention, then JDips answer already posted here fully applies. The resample and decimate commands provide the equivalent of anti-alias filtering with the additional good comment by Jdip on the questionable signal reduction by the high pass filter that the OP applied. Clearly we need more details on what the OP's true intentions are.
Perhaps not intended, but the OP has come up with a different approach to power spectral density (PSD) estimation that I haven't considered before (although it is quite similar to the idea of simply passing the FFT result through a FIR filter to smooth it). In interest of exploring that further, I will detail what the OP has done, how it relates and differs from existing power spectral density techniques (Bartlett and Welch methods); I will then provide some slight simplification to the OP's Python code, and then I will provide an example showing how the OP's results are similar to those existing methods on a test waveform using AWGN. In the end I have determined that the OP's PSD method has resulted in similar smoothing for the PSD but with more required processing to compute vs traditional methods (Bartlett and Welch's method). So the details provided here are just a fun study and further insight into power spectral density estimation but nothing more.
Fun Details (perhaps only for me)
The OP is not down-sampling in the time domain, but resampling the frequency domain samples themselves. Further we note that the OP has taken the magnitude of those frequency domain samples prior to using the built-in resample command to have less samples in the frequency domain. This doesn't change the sampling rate, but reduces the number of samples for a Fourier Transform at that same sampling rate.
If we for a moment ignore the magnitude computation, the operation done is effectively decreasing the time duration by a factor $D$ through applying anti-aliasing of possible time domain aliasing combined with selecting every $D$th sample in the frequency domain; this is the frequency domain equivalent to resampling in the time domain, which is effectively lowering the sampling rate by a factor $D$ through applying anti-aliasing of possible frequency domain aliasing combined with down-sampling by selecting every $D$th sample in the time domain.
This is very interesting, and not sure if the OP intended this, but by taking the magnitude first, the OP has provided yet another approach to computing the power spectral density with some similarity to the Bartlett Method that I haven't previously considered, in that the averaging of magnitudes results in less noise in the estimate of the power spectral density of the signal. This can be understood intuitively by understanding the following key points:
Samples of a real additive white Gaussian noise process (AWGN) in the time domain will be samples of a complex AWGN process in the frequency domain. By using Parseval's theorem and proper scaling in the DFT, we can see how the variance $\sigma_x^2$ of the real time domain samples can be equal to the variance $\sigma_X^2$ of the complex frequency domain samples. Further the variance $\sigma_X^2 = \sigma^2/2$ where $\sigma^2$ is the variance of the real or imaginary portion of the frequency domain samples.
The magnitude of a complex AWGN process is Rayleigh distributed; the mean of a Rayleigh distribution is $\sigma \sqrt{\pi/2}$, where $\sigma$ is the standard deviation of the real or imaginary portion of the underlying complex AWGN process. Thus, if we keep averaging the magnitudes of the DFT result (which is Rayleigh distributed, assuming a time domain AWGN), the result will converge to $\frac{1}{\sqrt{2}}\sigma_x \sqrt{\pi/2}$. Note that expressing this in dB results is the well known -1.05 dB under-estimate of the power when using the "square of the mean of the magnitudes" to estimate power of an AWGN signal instead of the "mean of the square of the magnitudes". This is known as
"true-rms" power detection for waveform agnostic power measurement (simple power measurement with electronics uses a diode to rectify and then a capacitor to average--- which provides the mean of the magnitudes, which is then sensitive to the actual statistics of the waveform).
What we do get by taking the magnitude first of shorter duration time domain samples, and then averaging is reducing the noise (the variance) of the the estimate for $\sigma$ (and the estimate for $\sigma$ is the average of the magnitudes as explained above). This is identical to reducing the Video Bandwidth (VBW) on Spectrum Analyzers (test equipment) which is simply filtering after magnitude detection, as opposed to Resolution Bandwidth (RBW) which is filtering prior to magnitude detection): Reducing the VBW reduces the variability in the noise floor from sample to sample (smoothing), but does not reduce the amount of noise as measured. JDip intends to update this post to explain more mathematically why the averaging of magnitudes results in less noise in the estimate of the power spectral density with the similar Welch Method.
What we lost by operating on shorter duration time domain samples is a larger resolution bandwidth which means less resolution in the frequency domain to distinguish between separate signals with similar frequencies. When estimating power spectral densities where the signal is sufficiently spread out over frequency (which is certainly the case for white noise), this disadvantage is of no consequence. We do care a lot about this when trying to estimate the power of very narrow band tones such as individual sinusoids or "spurs", as demonstrated in this post.
Until I find out that the OP's approach is already established, I will refer to the OP's approach as the "Tyesh Method" for comparison to the Bartlett and Welch's Methods of computing a power spectral density. A simulation result to be provided will also provide insights into any possible merits, if any, to using this approach.
OP's Python
As a minor comments regarding the Python used, the signal library uses Numpy which can operate on iterables (such as a list, Numpy array or other "list-like" collections) directly, so there is no need for the list comprehension. The following would provide an equivalent result.
downsampled = signal.resample(abs(fft(signal_chunk_filtered_new))[:int(chunk_size/2)], 1000)
However, in the interest of more closely modelling the processing in the Bartlett Method and to consider if it is a possible "low resource" approach to PSD estimation I modified this to be a simple moving average. Thus the resampling is simplified to be equivalent to a CIC filter, with the resulting "Tyesh Method" as I envision it to be a single FFT of the complete waveform, followed by a CIC decimator of the resulting FFT samples. (With proper gain scaling for resolution bandwidth and the -1.05 dB underestimate previously mentioned).
Is it really "low resource"??
No. This appears to just be a fun study and to confirm why Bartlett is still King. (and once we want to add windowing, for when we have a large enough dynamic range in the PSD to be concerned with spectral leakage, then Welch comes to the rescue as the "go-to").
With a quick study of required multipliers for each N pt DFT as $2N\log_2(N)$ we see that Bartlett comes out way ahead to the extent we can break a longer data set into shorter segments (at the expense of increasing the resolution bandwidth):
Bartlett: Qty N M-pt FFT's requires $2NM\log_2(M)$ real multiplies, then NM magnitude computations.
Tyesh: Qty 1 NM-pt FFT requires $2NM\log_2(NM)$ real multiplies, then NM magnitude computations.
Both require similar amount of additions, plus the additions for the post FFT averaging. Sure, I have come up with a way that allows Tyesh to do that averaging with super low resource CIC filters (and don't quite see how to pull that off yet with Bartlett, given the block by block averaging that is done in that case), but does that offset the increase in mults and adds in the FFT processing? I conclude note using the example demonstration below showing how this ends up:
Demonstration Comparing "Tyesh Method" to Bartlett Method
Below is the result of a simulation comparing the Bartlett Method for creating a power spectral density to the Tyesh Method. Note, this is using the same AWGN test waveform as demonstrated in the post also linked above, where we also see the result when using the Welch method and the significant noisy result we get when we simply take an FFT of the entire sequence with no further processing.
Visually, it appears that the Tyesh Method is a little noisier than the Bartlett Method, which surprised me as I really thought they would be similar: both are Rayleigh with the same distribution and averaged over the same number of samples, albeit different samples so will not match sample for sample but should have the same distribution in the averaged result. However, I measured the actual standard deviation of the power estimate (as plotted here in dB), and from that confirmed that the "noise" in both results are indeed nearly equivalent:
Standard Deviation of Tyesh Method Result: 8.15e-06
Standard Deviation of Bartlett Method Result: 8.03e-06
Here, no windowing was used with either approach, and for the Tyesh method the complete 100,000 sample waveform was FFT'd, and the resulting magnitudes were decimated by 256 with a first order CIC (this is equivalent to a 256 sample moving average followed by a down-sample by 256)...such a CIC is simply selecting every 256th sample from a single accumulator that rolls on overflow followed by a sample by sample difference at the output rate (so simple---- and importantly fixed point as RBJ would quickly point out if he was reading this.). For the Bartlett Method, the 100,000 samples were stacked into 256 point FFT blocks, each were FFT'd and then the magnitude of the results averaged block by block. Using the computation above for the total number of real multiplications needed just in the FFT processing alone for this example:
Tyesh: $2(100,000)\log_2(100,000) = 200,000(16.6) = 3,321,900$ mults
Bartlett: $2(100,000)\log_2(256) = 200,000(8) = 1,600,000$ mults
The Tyesh and Bartlett would each have the same number of additions as they have in multipliers in the FFT processing alone... I don't see the need to then review the savings in post averaging additions with the Tyesh method at this point and will just send it back into its corner of oblivion now that we've had our fun with it.
