What does "noise power is assumed to be the same at all nodes

Although noise is often called channel noise in various block diagrams of radio-frequency (RF) communications systems that communications analysts draw, the truth of the matter is that most actual RF channels are mostly noiseless, and what is referred to as the channel noise is primarily thermal noise in the resistive elements in the front end of the receiver.  The noise power in the receiver that communications engineers have to contend with is thus determined by the temperature of those resistors (and the value of those resistors).  Thus, assuming that the receivers in the various "nodes" are identical and also in the same geographical area so that all the receivers can be expected to be at the same temperature, it is reasonable to use a model in which the noise power is the same in all of the receivers.


Melodic or tonal pitch contour tracking?

I think a PLL would be hard to use for pitch detection in singing or speech. For the reasons you mention (classic phase/frequency comparators output garbage when the waveforms are too different due to a missing f0), but also because you'll need another method anyway for tracking the large changes in pitch or the switch to voice/unvoiced that would put a PLL off the tracks.

A baseline system could be to extract a pitch on small frames with a classic f0 detector (autocorrelation, AMDF, YIN, or any auditory model), and then use dynamic programming (or from another point of view, the Viterbi algorithm on a HMM whose states are discretized pitch values) to recover a pitch track. A ready to use Matlab implementation for this (auditory model + viterbi) is available here or (fancy autocorrelation flavour + dynamic programming) here.

This will be discretized on both the time and pitch axes, but you could then interpolate.


Ramp function as derivative in frequency domain?

You may be familiar with the classical properties of the Fourier transform. One of them is the one regarding differentiation in the time domain:

$$x(t)\xrightarrow{\mathscr{F}}X(j\omega) \implies \frac{d}{dt}x(t)\xrightarrow{\mathscr{F}}j\omega X(j\omega)$$

And there appears the ramp function, $j\omega$.

The demonstration of why this is true is straightforward. If we have a signal $x(t)$ with Fourier transform $X(j\omega)$:

$$x(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}X(j\omega)e^{j\omega t} \ \mathrm{d}\omega$$

Then by differentiating both sides:

$$\frac{d}{dt}x(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}j\omega X(j\omega)e^{j\omega t} \ \mathrm{d}\omega$$

And that's the (pretty simple) proof.


Assuming that $N$ is the length of your signal $s$, the normalized signal $s_n$ is given by:

$$s_n = \dfrac{s}{\sqrt{\dfrac{\sum_{i=1}^{N}\left|s_i^2\right|}{N}}} $$

The denominator is nothing else than Root Mean Square value of your signal. Thus the code is doing a simple RMS normalization.

You can think of it as a method of normalizing the average of your signal and still allowing for some peaks to clip it (instead of being set to 1). In other words the standard division by the maximum absolute value of your signal will always guarantee that the sample values will be within the range of $[-1;1]$, whereas RMS normalization doesn't.

This method is widely used for audio processing and speech processing.


Does a matched filter really improve SNR?

There is no contradiction. Matched filters are often used because they are an implementation of a sliding correlator, which is the optimum detector (in the maximum-likelihood sense) for linear modulation in the AWGN channel. That is, the SNR at the output of the matched filter will be maximized when the filter lines up with the point of reception of the signal that the filter matches.

This doesn't contradict the known properties of your receiver (i.e. its noise figure). The point to notice is that the SNR can be different at the input of the matched filter than it is at the filter output. The SNR at the matched filter input (i.e. your front end's output) is determined by the signal power at the receiver input and the noise figure of the front end. But by applying signal processing techniques, like matched filtering, you can improve that SNR to achieve more accurate detection.


Matlab FIR filter numerator coefficients to python

You don't want np.ones_like(num) in the denominator.  That's making a denominator of all ones, with lots of possibly unstable poles.  Your denominator in Python needs to be [1,0,0,0,0,0,...]. Or just [1] if possible.


The easiest proof uses the convolution property of discrete LTI systems.
for Continuous time:
If you write down your periodic sigal as it's Fourier series (with omega =2pi/period) and describe the output using the convolution property you have:

We get the output as Fourier series with the same period.
now the exact same thing can be applied for discrete system (by using the discrete convolution).

2)What conditions are required for this system to be linear?.

difference equations are analogous to differential equations. With the initial conditions at rest (if x(t)=0 for negative t then y(t)=0 as well) The differential equation describes LTI system.
The same can be shown for the discrete time case.
With initial rest Your equation will describe LTI system.


Why don't DSLR cameras default to png instead of jpeg?

Its force of habit and good enough. People were using JPEG's in their digital cameras before, so they stuck with it, and its plenty fine for viewing (and it produces pretty small file sizes as well). You're free to turn to some other compression method your camera supports.

Note that the RAW files often contain a lot more information than what a PNG would have anyway (like a bunch of stuff about the sensor), which can be exploited in programs like photoshop for editing, so from that perspective, you may as well have the RAW rather than the PNG (especially since the RAW is rather space efficient to begin with). You can also put in things specific to your camera there if you want.

Note that many cameras have the option to shoot RAW+JPEG so you simultaneously get both, which is good for viewing (use the JPEG) and editing (use the RAW).


If I understand the question properly, then it sounds like the difference between time-domain multiplication and frequency domain multiplication.

Windowing is a time-domain multiplication operation.  As such, it is equivalent to convolving the spectrum with the window kernel, which will operate all frequencies the same way.  On the other hand, filtering is (can be) a frequency domain multiplication, so the filter can be designed to operate specifically on chosen frequencies (or bands of frequencies if you like).

To look at it slightly differently: if you consider your signal as the sum of all your subcarriers (which you can do with linearity), then windowing your signal with an RC is the same as windowing each of your subcarriers with an RC before adding them together.


Conversion of airflow signal into 1-D as a input for CNN

The airflow signal is usually respiratory input in a polysonogram. If your input files are .EDF files, then you need to extract the airflow signal only, which is a 1-D signal. I think there are many implementations of edfread.m that are available online.

https://www.mathworks.com/matlabcentral/fileexchange/31900-edfread


Isn't Gardner's algorithm and Early-Late gate the same thing?

The Gardner Timing Error Detector is diagrammed in the graphic below, where two samples per symbol are used, and the error is determined using Prompt*(Late-Early), and when synchronized the center sample (Prompt) will be midway between two symbols. In contrast, an Early-Late approach uses (Late-Early), typically on a correlated symbol response, and when synchronized the Prompt sample which is between Early and Late will be centered on the correlation peak.

Note that the Gardner Timing Error Detector forms a Timing Error Discriminator by approximating the derivative of the waveform at the assumed zero crossing and weighing that by the actual sample value at the crossing location. (The difference equation $y[2n+1]-y[2n-1]$ is an approximation of the derivative).

In contrast, a Early-Late gate approach approximates the derivative at the peak correlation of the data symbol. When timing is synchronized in a Early-Late gate approach, the two Early and Late samples will be a fraction of the data symbol on each side of the peak correlation. They could be spaced 1/2 a symbol, making the spacing the same as the Gardner, but the key distinguishing feature is the lock point of the center sample.

The functional operation of the Gardner as forming a discriminator by multiplying the waveform with its derivative is demonstrated by the two graphics below.

The functional operation of an Early-Late approach is shown below for the case of timing being synchronized, too late and too early. A GPS receiver is a good example where this approach is often used, where the output of the correlator of the received sequence with the GPS acquisition code used for every symbol would have the triangular response shown below, which specifically is the magnitude of the correlator versus timing offset with the code. Three correlators would be used, with the Early and Late correlator spaced approximately one half of the correlation width as shown in the diagram. The Early and Late correlators would be used to keep the Prompt correlator centered on the peak, which would be used as the recovered symbol.

Notice the key difference between the two approaches. In the Garder Loop the derivative is fairly constant in vicinity of the lock point, and it is the waveform itself providing the discrimination curve (change of error magnitude versus change in timing offset) - while in contrast with the Early-Late approach it is the derivative itself that provides the discrimination curve.


How to calculate the impulse response from the transmission function of an optical device?

If the function $t(x,y)$ is an optical transfer function (OTF) then your goal is not to find a function $h(x,y)$ such that the convolution with your input $f(x,y)$ equals the multiplication with $t(x,y)$. Convolution with the impulse response (point spread function (PSD) in optics) happens in a different domain. The OTF is the Fourier transform of the PSD, and the equality

$$f'(x,y)=f(x,y)t(x,y)$$

corresponds to the equality

$$\mathcal{F}^{-1}\{f'(x,y)\}=\mathcal{F}^{-1}\{f(x,y)\}\star h(u,v)$$

where $h(u,v)$ is the impulse response (PSD) given by

$$h(u,v)=\mathcal{F}^{-1}\{t(x,y)\}$$


Sparse CWT matrix-shape data lossless compression

Start simple: that calls for any Run-Length coding (RLE).
Also, floating point makes very little sense for image data, so the first step that will probably save you 50% of space is conversion to unsigned fixed point (multiply every value with 2¹⁶, then divide through the largest value of the original image, round to nearest integer, keep the result as uint16).
5% mean absolute error means: you don't need floating point; rounding small values to zero is absolutely¹ acceptable to you. So, drop the float idea and convert to integer, as described above.
If for some (somewhat untypical for image data) reason you want to keep it floating point, convert it to a format where you store the mantissa and exponent bits separately and the latter at all only for non-zero mantissa values. Since pixel data typically only has a positive sign, you don't even need to store the sign bit.
After reducing the number of exponents you've stored, maybe any further compression isn't necessary to achieve your goal (don't overengineer!).
Then, look at your data: It's also quite likely that your exponents don't need the full 2⁻¹²⁶ to 2¹²⁷ range, so dropping bits on your exponents makes a lot of sense, too. Make sure you're not totally killing it with the number of mantissa bits, too. Wild guess: half as many as in 32 bit IEEE754 will do.
If after doing these thing it still is necessary to compress further, RLE or Lempel-Ziv would be appropriate solutions.
(LZ4 is kind of exotic; the more established LZMA algorithm as e.g. in XZ should work sufficiently well, and it really doesn't sound like you're CPU bound here...)
In all seriousness, though, sounds like a job PNG, high-quality JPEG, or if you have access to that, JPEG-XL, could do out of the box, so I'd simply look into converting my data to a format that makes sense to readily available lossless image codecs instead of trying to roll my own.

¹ hurr durr, a pun!


Deriving statistics of band limited Random Noise

at first read a book on stochastic process or stochastic signal processing.
it's better that we say we have a random process which it's first order distribution (PDF of the process at each instant of time) is independent of time and is a Gaussian with known mean and variance. Now we want to obtain the first order properties of filtered process.

To begin we write the output of filter as a convolution. The output process at each instant of time is a linear combination (weighted average) of input process at different instant of time which all has the same Gaussian distribution, also we know linear combination of variables with Gaussian PDF leads to another variable with Gaussian distribution. So the first order distribution of output must be a Gaussian.

Now taking expectation from both side of convolution, we see expectation of filtered process at each instant of time becomes convolution of filter response and the mean of input process which is a constant over time, so it becomes the multiplication of input mean and the DC response of filter.
To obtain the variance of process around its mean, we have to find the expectation of squared difference of the process and its mean or expectation of squared value of output minus the square of it's mean. To relate the squared value of output to the input we could convolve the filterd signal with its reversed time version at zero lag. At the end we see new variance is the variance of input times square of DC response of filter.


confusion between frequency and amplitude in digital to analog and vice versa conversion

The discussion of harmonic in this context is to show with a simplest approach the maximum bandwidth needed to contain a digital bit stream.

The highest frequency for a digital bit stream occurs when we transmit a sequence that changes on every bit:

1 0 1 0 1 0 ...

with the analog representation of this being a square wave. So if we were to use a simplest D/A converter, the analog output in this case would be a square wave IF we had infinite bandwidth in the frequency domain.

What is useful at this point is to draw on the Fourier Series Expansion, where we have learned that a square wave can be decomposed into an infinite sum of frequencies, starting with a base frequency which is a sine wave that exactly fits into the 0 1 0 1 0 1 0 pattern, which we call the "first harmonic". For a general Fourier Series Expansion, the theory continues that all continuous periodic functions are decomposed into such a first harmonic that is at frequency 1/T where T is the time duration of the base function (in this case 0 1 as a square wave) AND all the higher integer multiples of this frequency (the 2nd, 3rd, 4th.. etc harmonics) some of which may be zero, but if the pattern is repeating (periodic) there will be no frequency components anywhere else. In the case of a square wave it is specifically all the odd harmonics 1, 3, 5... that if you add up more and more of these (each with an amplitude that goes down at 1/f or 1, 1/3, 1/5, etc) the result will get closer and closer to a square wave. You need very high frequencies to represent sharp (fast) edges in a time domain signal!

For the purpose of transmitting information, the pattern 0 1 0 1 0 1 CAN be transmitted with a simple sine wave, as we can use a threshold detector at the other end to reliably convert it back to the 0 1 0 1 pattern and thus recover the information.  As we increase bandwidth to include more of the higher harmonics, the signal will get closer and closer to a square wave.

So they are just making the point here that you can use this logic to determine what is the minimum bandwidth required for a digital bit stream. 0 1 0 1 is the highest frequency pattern and we showed above that we can send this with just the first harmonic. Any other patterns such as 0 0 1 1 0 0 1 1 0 0  contain lower frequency components; so as long as our maximum bandwidth exceeds this, these can be transmitted without error as well. In the end our random data stream can be viewed as multiple patterns that are summed together, so in frequency instead of seeing discrete tones, the energy will be spread evenly across that bandwidth (if our data is truly random with one bit completely independent of the next).

I hope this helped!


DFT vs FFT: does odd number works for DFT only?

An FFT is an optimization for computing the DFT that works well for any number of samples that is the product of mostly small prime number factors.  They produce the exact same results, except for numerical noise. The factor 2 is common, but 3 and 5 are also small primes that can be used in a very efficient odd number length FFT implementation.  Many newer FFT libraries (Accelerate/vDSP, FFTW, etc.) thus support odd length FFTs, not just powers of 2.  On newer processors, the speed of an FFT is more limited by data cache size and policy and CPU instruction pipeline depth and hazards, not by the input being only a power-or-2 or not.

The algorithm makes no assumption about the signal being periodic (although some people do), as it is just a finite matrix basis transform.  If the signal input isn't exactly integer periodic in the FFT width, with that width being cut out of a longer signal, then there will be some windowing artifacts, often called spectral leakage (but they have the exact shape as the transform of the window used, so they are far more related to the window than to the spectrum.)

If you change the size of a DFT or FFT window so that it is an exact integer multiple of the period of a purely periodic input waveform, then the windowing artifacts (or leakage) just happens to be all zeros.


how does the SDR sample rate translate to sample rate and bit depth of the resulting audio signal?

Extremely indirectly. I haven't got precise mathematical formulas for you, but here's an overview of the topic which should let you figure out where you want to dig deeper.

First of all, FM is an analog modulation. This means that there is no inherent sample rate or bit depth of the signal recovered by demodulating, but only an audio bandwidth and noise floor. So let's look at what influences that:

The bits per sample (bit depth) of the digitized RF determines the amount of error introduced by the quantization (conversion to discrete samples), which can be described as quantization noise.

Note that quantization noise is at a constant level in the digital signal coming from the ADC. This means that — unlike radio interference — the signal-to-quantization-noise ratio improves the more you turn up your receiver's gain. But you can only do that so far before you get clipping (signal level exceeding the representable sample values).

Therefore, the best results are obtained when you have a narrow analog filter in your receiver to remove all but the desired signal before digitization — the wider bandwidth you receive, the more total power is arriving, and so the less gain you can use before clipping.

Of course, using a narrow bandwidth prevents you from having those pretty waterfall displays of ten stations at once. It's a tradeoff.

For more words on the subject, see my previous answer over on Amateur Radio Stack Exchange, to the question Why is dynamic range relevant for an SDR?
A simple analysis would suggest that the sample rate of the digitized RF has no effect, because it's just bringing in a wider bandwidth, full of unrelated signals and noise, than we actually need. However, this is not true, if the ADC has fewer bits per sample than the actual DSP processing (for example, 8-bit integer samples from the ADC going into 32-bit floating-point DSP operations). The extra samples contain some of the information that was lost in quantization; resampling the unnecessarily-high-rate signal to a lower rate will automatically use that information, so those 32 bits aren't actually wasted.

So to combine this with the previous point, the best case given fixed hardware is to configure it (if possible) for the narrowest analog filter bandwidth and the highest sample rate (creating a signal whose spectrum has an obvious filter hump in the middle) and then resampling it before demodulation.

To sum up, the sample rate and bit depth of the incoming RF determines how much noise is added to the signal by the process of digitization. The next question is how that noise, the RF noise entering the software demodulator, affects the noise level of the resulting audio, and I'm afraid I don't know that and haven't found an obvious reference.

Now that we know that our input can be considered as a signal with more or less noise depending on its parameters, let's look at what to do with the output.

As Nyquist tells us, the sample rate of the output audio signal does not need to be greater than twice the highest frequency present in the signal. That's the audio bandwidth, which for broadcast FM is, let's say, 15 kHz. Therefore, the required sample rate is 30 kHz, ideally. However, the signal must be filtered to remove noise that is outside the expected bandwidth, and using exactly 30 kHz sample rate means that the noise which passes the necessarily-imperfect filter will be aliased down into lower frequencies. Thus, it is better to use a higher sample rate to allow some room for the filter transition band — perhaps 44.1 or 48 kHz, because those are common audio sample rates.
Quantization introduces quantization noise — but we already have noise in our signal. Therefore, the bit depth of the audio should be chosen so that it introduces an insignificant amount of quantization noise compared to the noise already present. (Note: You can't avoid introducing any noise at all. Your finite-precision arithmetic is doing that at every step.)


FFT calculates spectrum over the whole duration of your signal (integrates energy at each frequency over the time of accumulation) while spectrogram shows frequency variation in time (instantaneous frequency). Divide the whole frequency band of interest in some equal intervals and evaluate how much time signal's frequency falls in each interval. These time periods aren't equal. So by FFT frequencies corresponding to larger time periods gains more power and will be major and vice versa. Now lets look at your spectrogram. We see that it's about 100 and 200 Hz when the most stationarity in frequency exists, there are maximums of the spectrum. Frequencies between 100 and 200 Hz appears 2 times, so we see raising of the spectrum at this span. Low and high frequencies appears only 1 time and frequency at these regions varies faster: the spectrum after FFT decreases.


Yes, it is possible to determine the path that the Viterbi algorithm found through the trellis. Just apply the encoding algorithm to what you are calling Y_dec and you will get the corresponding codeword of length $2048$. You can then compare it to Y, the transmitted codeword, to see where the channel made errors.
Additional notes:
If the data to be transmitted 1024 information bits and the code to be used is a half-rate $(2,1)$ convolutional code with constraint length $7$, then the convolutional codeword produced by the convolutional encoder is $2060 = 2048+12 = (2\times 1024) + 2\times 6$ bits long where the extra $12$ bits are produced by the encoder as zeroes are pumped into the $6$-bit-long data buffer as the buffer empties out. I don't speak MATLABian gobbledygook at all and so don't know what Y_dec(1+tblen:end) means, or indeed why tblen is involved at all.  If  Y = convenc(X, trellis) produces a codeword vector Y of 2060 bits, then
Y_dec = vitdec(Y, trellis, tblen, 'cont', 'hard') should produce X again. I don't see where there was any noise involved and so the Viterbi decoder should just zoom through the trellis following the same path as the encoder so that Y_dec should equal X. If there is noise, then Y_dec should be the same as X in most of the 1024 bits except for some bursts of errors in that span. In this case, convenc(Y_dec, trellis) will not be the same as Y, but in the absence of any channel errors, Y_dec=X and  convenc(Y_dec, trellis)=convenc(X, trellis)=Y.


Matched filter for recieved RF digital packet by LPF

Well, you know that the output of your filters shouldn't change faster than a symbol period. You could low-pass filter these to a bandwidth only a bit bigger than the symbol rate.
Since a filter is a linear operation, you can commute it and have it before the matched filters, so you only need one low-pass filter.

smoothly and to do zero crossings in order to choose when I have 0 and 1 in binary

Don't. You should be using a threshold instead of going by zero crossings: the zero is your maximum likelihood decision boundary.
You'll need symbol timing recovery, anyways. Since this is a bursty system, not an "endless" stream of symbols, I'd guess the packet preamble is probably meant to give you that information – use it!
If you can correctly estimate symbol timing from the preamble, your Matched Filters can be coherently decimating down to the symbol rate.


Sampling at a higher rate will distribute the quantization noise over a wider frequency, thus reducing the noise spectral density due to that quantization noise, with a lot of caveats. For more details on that see
What are advantages of having higher sampling rate of a signal?

In your last paragraph, if you are referring to running the same filter at a lower rate versus a higher rate, at a lower sampling rate the digital filter will require LESS taps for the same filter performance. For more details on that see Filter Order Rule of Thumb and How many taps does an FIR filter need?.

An effective strategy to consider is to use a higher sampling rate which relaxes the analog filtering requirements and can give you more effective bits if needed at the analog to digital boundary (which are realized after subsequent filtering digitally). This is then followed by efficient resampling techniques in the digital domain to get to a lower sampling frequency prior to providing the final filtering as required in your system design (coarse and simpler filters are done at the higher rate, and then final "shaping filters" are done at the lowest rate possible thus minimizing number of taps, power dissipation and resources).

For a further details on decimating to a lower rate for improved filter performance see Fast Integer 8 Hz 2nd Order LP for Microcontroller.


All depth measurements I've been able to read are very accurate (±0.1cm)

That's highly unlikely. It seems that something is quantizing your data at pretty rough interval: maybe 2cm or so plus some extra noise. There are is also a spike at t = 2000 units ?) where the depth jumps from 0.1m to about 0.5m in one sample. Your sample interval is 40 milliseconds. 0.5m in 40ms corresponds to a speed of 12.5 m/s. Can your object really move that fast ?
Overall your data looks quite bad (noisy, quantized & spikey) and not like a real physical object tends to behave.

find out the speed of the body by differentiating the measured depth value

The second curve doesn't look the derivative of the first curve at all. I would expect a massive velocity peak of 12.5 m/s at t = 2000, but your velocity curve just shows a minor blip. Another big peak should be at t= 4000. I have no idea what you did there, but it doesn't appear to be differentiation.

This is a simple calculation of the "velocity" from the depth measurements using Excel. It seems more reasonable than the plot in the question.
Note that the data used was from Web Page Digitizer, so it won't be quite the same as the original.

In the future please label the X and Y axis of your graph with proper units. I have no idea what the unit of your velocity are supposed to be but it's unlikely to be in m/s.


Design of the low pass filter when decimating

You can determine the "best" filter by analyzing what type of errors is created and how much of each error your specific application can live with.

If your original content is already perfectly band limited, than you don't need a low pass filter at all and can directly throw away the extra samples.

The decimation low pass filter typically determines how the band limiting of the content is done. Here are the factors to consider

Aliasing: since the low pass filter is never perfect, the decimation will
create some amount of aliasing. The amount depends on the steepness of the filter and the energy content of the signal at the frequencies that need removing.
Pass band: The low pass filter isn't perfectly flat in the pass band so there is some amount of ripple you need to live with. The pass band must be below the new Nyquist frequency. The higher the filter order, the closer you can get and the wider the pass band will be.
Phase response: Any IIR filter will change the phase response of your signal. The steeper the filter, the higher the maximum group delay. If need to maintain the phase, you need to use a linear phase FIR filter (ideally in a polyphase configuration).

The choice of filter is a trade off between pass band width, pass band amplitude distortion, pass band phase distortion, aliasing noise and computational complexity.If you have a spec for all of these, you can design the "best" filter. An example would be "pass band +-0.5dB up to 90% Nyquist, aliasing at -80dB or below, linear phase".


Smoothing in frequency domain and impulse reponse of the filter in time-domain

You get the time domain values of hamming window from the hamming function. Multiplying it with the frequency domain data will not give the correct result. Smoothing is usually performed by multiplying in time domain. You could try multiplying the fourier transform of hamming window with the Frequency domain signal. and then performing the inverse fft.


Intuitive interpretation of Laplace transform

Why is the fourier transform a special case of the laplace transform?

The Laplace transform produces a 2D surface of complex values, while the Fourier transform produces a 1D line of complex values.  The Fourier transform is what you get when you slice the Laplace transform along the jω axis.  For instance, a simple lowpass filter $H(s)=\frac{1}{s+1}$ has a single pole in the S plane to the left of the origin:

Viewed from the side, the magnitude of this Laplace transform forms a surface, with the pole acting like a tent pole that raises the amplitude to infinity at that point (and an implied zero at infinity that drops the amplitude to zero the farther away from the origin you get in any direction):

If you now take the value of the surface along the jω axis only, that's the Fourier transform.  It's the red curve in the image above, which you can see forms a lowpass filter.  If you moved the pole farther away from the origin, the tent would move in the same direction, and the slice along the jω axis would drop, both reducing the gain (which we compensate for by adding overall gain) and increasing the cutoff frequency.  I've been meaning to make some animations of stuff like this...

http://www.maximintegrated.com/en/app-notes/index.mvp/id/733

https://dsp.stackexchange.com/a/9579/29


How to represent multiple related signals via mathematical function?

Try the Eureqa software. It tries to find the simplest mathematical relation between multiple data. If the relation truly has a simple symbolic form, I think it should find it. Otherwise it will produce some strangely formulated approximation.


Are there other filter models besides IIR and FIR for processing digital signals?

FIR and IIR filters are discrete-time linear time-invariant (LTI) filters. There are also (linear) time-varying systems, which could be implemented just like FIR or IIR filters, but with time-varying coefficients. Audio effects like chorus, flanger, or phaser are time-varying filters.

Of course there are also non-linear filters, the most well-known of which is probably the median filter, which is very useful for removing impulsive noise. Many different types of non-linear filters are used in image processing.

Note that none of the above mentioned filters can be characterized by their impulse response or, equivalently, by their frequency response. This is only possible for LTI systems.


Is it required to quantize wavelet coefficients before one derives features from?

The global answer is no, it is not mandatory. Yet weak forms of quantization (of values and locations) may help, like absolute values and range compression.

First, wavelet coefficient values can be complex. Wavelet  transform are usually linear, and are often shift-variant.  If you want a classification that does not depend on data sign or phase, and on (some) offsets, ie if you want the same class for:

black object on white or white object on black.
slightly offset signals

at least you need to take an absolute value, one of the weakest form of quantization.

Then, depending on what classification you are used, it is not mandatory. You can for instance fit a distribution to wavelet coefficients in each subband, and use the distribution parameters for clustering, without further quantization. For instance, generalized Gaussian-Laplacian distributions were once considered a good fit for wavelet coefficients, see for instance Wavelet-Based Texture Retrieval Using Generalized Gaussian Density and Kullback–Leibler Distance. The following plot is taken from a paper related to adaptive filtering, however this generalized Gaussian modeling from wavelet frame decomposition with different power laws may apply to classification too.

More generally, range compression, like applying a power-law (Box-Cox, Anscombe transformations) to absolute wavelet coefficients, is often used to set thresholds, that are another form of classification.

The last quantization we often do not think about is location. We previously talked about values only. But you can gather wavelet coefficients in vectors or blocks, packed into a "Reference Observation Vector", and compute statistics on them, that will be used for further processing.


For an undergrad DSP lab, there's absolutely no benefit in complicating what hardware your students need to learn: Let them do DSP on a completely normal PC.
There's plenty of offline processing and real-time capable DSP platforms for PCs. I'd personally recommend getting a bunch of Linux PCs, and having Python3 with scipy/numpy/matplotlib on them for basic algorithms exploration on that.
For advanced topics, like doing live signal processing on audio coming from a microphone or from an internet stream, or doing software defined radio with real SDR hardware, I'd go for GNU Radio, which comes with a lot of useful building blocks, and can be easily extended using Python and C++.¹

¹ Disclaimer: I'm involved with the GNU Radio project, but I do, and I know several educators, use it in DSP/communications engineering labs. Say "hello" on the mailing list.


Concept About Estimated Standard Deviation

Given data $ { \left\{ {x}_{i} \right\} }_{i = 1}^{N} $ the Empirical STD of the data is well defined:

$$ STD = \sqrt{ \frac{1}{N - 1} \sum_{i = 1}^{N} { \left( {x}_{i} - \bar{x} \right) }^{2} } $$

Where $ \bar{x} $ is the empirical mean of the data given by:

$$ \bar{x} = \frac{1}{N} \sum_{i = 1}^{N} {x}_{i} $$

Now, if there's a model on the data (Such as signal and noise with certain CDF) the empirical calculation should be altered accordingly.

For instance, given a signal which is linear with AWGN the STD of the noise can be estimated by removing the linear estimated signal first.


The first and best solution is to increase the input gain which should be possible on most platforms. If not, then read on...

If multiplying by a constant caused distortion then either:

it is too loud and distortion is inevitable (barring what I'll just call "fancy" signal processing)
you are doing something wrong.

I'll assume #2 since you started with the case of "barely audible", and you should learn how to apply gain before doing anything "fancy" anyway. To apply gain, make sure you multiply each sample individually by the same constant. Remember that the sample count = frames * channels. Loop through each each sample in the buffer and multiply it by the gain.

For this case, you probably want to use the float datatype in PortAudio, but any type will work.


Python toolboxes for state-space estimation via subspace estimation

You can use pyvib to do frequency based subspace identification. Beware that there is no estimation of the initial state. It is possible to do optimization of the identified model, if the data is not perfectly linear. See the implemenentation, maybe you can use it, in case you want to do your own implementation.

Somewhat incomplete example. Take a look at the example to get a working code:

git clone --depth 1 https://github.com/pawsen/pyvib.git
export PYTHONPATH=pyvib

from pyvib.subspace import Subspace
from pyvib.signal import Signal
# partion the data so the format is (npp,inputs/outputs,R,P)
sig = Signal(uest,yest,fs=fs)
sig.lines = lines  # which freq lines should be used
sig.bla()  # computer best linear approx.
nvec = [2,3]  # model size to scan over
maxr = 5  # max number of rows
model = Subspace(sig)
models, infodict = model.scan(nvec, maxr)
errvec = model.extract_model(yval, uval)  # extract best model on fresh data
model.estimate(n, r)  # or do direct estimation if you know sys size

I have written the code. You are more than welcome to criticize and suggest improvements.


Carrier sense multiple access, MAC layer

In collistion detection, a node may start transmitting at any time, and backs off if it detects that it caused a collision. This forces the node that was already transmitting to back off, too. Its advantage is simplicity; its disadvantage is that it introduces big inefficiences: a packet makes it through the network only if no other node decides to transmit at the same time.

In collision avoidance, a node first looks at the medium. If it detects a signal, then it backs off and waits before attempting to transmit again. This increases the transmitter's complexity a bit, but the network efficiency is increased. A collision only happens when two nodes decide to transmit at (almost) the same time.

You can find the details and how to calculate the specifics in any textbook on networking.


Similar to a a vector that is normal when it's magnitude is 1, which is the inner product with itself. In the same fashion, a function is normal over a defined range $[t_1,t_2]$ when the root of the integrated squared magnitude is one (if it's normal the root is trivial but this shows the same form as derived from an inner product):
$$||x(t)|| = \sqrt{\int_{t_1}^{t_2}|x(t)|^2dt} = 1 $$
Normalization of functions is done with a simple scaling: If the above formula results in a number of than one, then the entire function is scaled by that number.


Why do the two methods give different answers for the Fourier transform of $Y = \cos(\omega_0 t + \phi)$?

They are the same.

Using the Dirac impulse sifting property $f(x) \delta(x-a) = f(a)\delta(x-a)$ you can also verify that the second method produces the same output as the first method.


Determining the Nyquist rate of signals added and multiplied together

Let $Y(\omega)=X(\omega)+X(\omega)e^{-j\omega}$ with $X(\omega)$satisfying $X(\omega)=0$ for $|\omega|>\omega_m$. Then, by its definition, $Y(\omega)$ also satisfies $Y(\omega)=0$ for $|\omega|>\omega_m$, i.e., the maximum frequency remains unchanged, and, consequently, also the Nyquist frequency remains unchanged.

Note that if $x(t)$ is the signal corresponding to the spectrum $X(\omega)$, then $X(\omega)e^{-j\omega}$ corresponds to the signal $x(t-1)$, which clearly has the same maximum frequency as $x(t)$.


Let's assume we have a sample rate of $f_s=10 kHz$ and FFT size of $N=1000$

Your bin spacing is $\delta f = 10 Hz$. It's simply the sample rate divided by the FFT size. That's all there is to it.
Keep in mind that the FFT requires both the time domain and the frequency domain signal to be perodic. Most signals are not, so there needs to be some sort "accommodation" that's specific to the needs for your application.
A real-valued sine wave with a frequency that is an integer multiple of the bin spacing will only show up in two bins symmetric around zero. I.e. a 100 Hz sine wave will show up in bins +10 and -10 (assuming we index the bins from $[-N/2,N/2-1]$
A complex exponential  with a frequency that is  integer multiple of the bin spacing will only show up in one bin. I.e. a +100 Hz complex exponential will show up bin +10.
A cosine wave with a frequency that's NOT an integer multiple of the bin spacing (e.g. 105Hz) will spread out over all bins. The bulk of the energy will be in the neighborhood of 105 Hz but there will be energy in ALL bins. This is a direct consequence of bullet 2): the FFT assumes periodicity with $N/f_s$. In essence the FFT doesn't calculate the spectrum of the sine wave, but the spectrum of a signal that consists of 1000 sample long snippet of a 105Hz cosine that gets repeated infinite times. At the splice point there is a massive discontinuity and that causes the energy to spread over all bands
In order to minimize this spreading effect, one can apply a window, that tapers down the edges of the time domain frame and reduces the discontinuity. On the other hand, the window also smears out the energy in the frequency domain, so you can't get any sharp lines any more.
Many people think that the FFT will give you an "exact" representation of the signal's frequency content. That's not the case unless the time domain signal is either periodic with the FFT length or if it's finite in time. For all other signals the results will depend heavily on the exact process and parameters chosen.

EDIT: reworded item 7 per the comment


Deriving the Convolution Kernel of the Inverse of a Signal

Since Convolution is a Linear Operation the easiest way to derive the coefficients is using a Matrix Representation of the problem.

Take the kernel $ y $ and build from it the Convolution Matrix.
Then write the vector representation of the Delta Funciton $ d $ and solve:

$$ Y x = d $$

In case the problem is solvable (Namely there is an inverse operator of $ Y $) you'll get you answer.
Of course you can always use Least Squares Solution.

Full code for building Convolution Matrix was given in Estimating a Signal Given a Noisy Measurement of the Signal and Its Derivative (Denoising) and Deconvolution of Synthetic 1D Signals - How To?


The magnitude of the roots does not give you the format gain. You need to evaluate the frequency response of the filter $H(z)=\frac{1}{A(z)}$, where $A(z)$ is the polynomial returned by the LPC function, at the point on the unit circle corresponding to the desired format frequency. Like this freqz(1,a, angle(r)) for instance.


I believe this will do:

x=a.^((0:40)+(0:40)’)


Playing/Recording acoustic signal in matlab to tell difference between 2 sounds

I understand that you want to tell if your second laptop is placed on the left or right (or maybe on the axis) of main node?

I suggest you then to play the two-tone stereo signal, i.e.: $15\;kHz$ from the left, and $16\;kHz$ from the right speaker at the same time. Assuming that you are going to use sampling frequency of $f_s=41000 \;  \text{Hz}$ and taking the average distance between two speakers in a laptop to be $25\;cm$, you can calculate that the maximum delay between two signals, if the second node is on the right, will be approximately: $\dfrac{0.25}{343}=0.7289\;msec$. This corresponds to $32$ samples at a given sampling frequency. It is therefore a hard task, but still possible. Then after finding the delay between two signals, from its sign you are able to tell which speaker is further away, ergo you are on the opposite side.

In MATLAB you just need to generate proper signals, given sampling frequency and play them back/record with: MATLAB - Record and Play Audio.

Few problems can arise:

You want to use very high frequencies, and this can be problematic as laptop speakers are usually low-grade gear and they possibly might have some problems- same goes to microphone. Also remember to turn any hardware noise cancellation, spatial enhancements, etc.
Detection of you signal onset could be problematic as the resolution of DFT will be to low. You can also use some kind of correlation detection, but this is a problem for later on.


Equation that describes a gravity drained tank system

Gravity-fed outflow is proportional to the square root of the liquid level (see Bernoulli's equation ). Let's call $h_i$ the height of the upper (inflow) tank, and $h_o$ the height of the lower (outflow) tank. $q_{12} = k_{q_{12}} \sqrt(h_i) $ and $(dh_i/dt) = k_{h_i} (q_i - q_{12})$.

For the system outflow, $q_o = k_{q_o} \sqrt(h_o)$, and $(dh_o/dt) = k_{h_o} (q_{12} - q_o)$.

Valve assumptions: linear and fast. So the steady-state flowrate of $q_i$ is proportional to the valve opening via $k_{v_i}$, and the dynamics of the valve are much much much faster than the tank level dynamics, so no differential equation is needed - algebra only.

Under P-only control, $q_i = k_{v_i} k_c (h_{o_sp} - h_o)$

You might want to linearize those square root terms via Taylor expansion to make the Laplace transforms a bit less hairy, if you're working in the $s$ or $z^{-1}$ domain.

Of course, in industry (and assuming that the level of this tank is not critical - it rarely is) you'd just set it up to P-only control with a dimensionless gain $k_c = 1.67$. This way, that if the setpoint is 50%, the valve is wide open when the lower tank level is at 20%, and fully shut when the lower tank is at 80%. And then you'd move on to the next crisis.  :)


What I was getting at in the comments above is that the linear system $w_1(t)a(t) + w_2(t)b(t) + w_3(t)c(t) = k$ has an infinite number of solutions, so you need to state some criterion that allows you to choose a unique solution. I think you have pointed out a constraint that is worth examining.

The idea: note that the linear equation I gave above is the equation for a plane in $\mathbb{R}^3$:

$$w_1(t)a(t) + w_2(t)b(t) + w_3(t)c(t) = k$$

This defines a plane in three-dimensional Euclidean space.
$\left[a(t), b(t), c(t) \right]$ is a vector that is normal to the plane.
The collection of all vectors in $\mathbb{R}^3$ of the form $\left[w_1(t),\ w_2(t),\ w_3(t)\right]$that satisfy the above equation are points on the plane.
Since you said that $w_1(t), w_2(t), w_3(t)$ don't change much over a short time duration, then you can assume that the point $\left[w_1(t+\Delta t),\ w_2(t+\Delta t),\ w_3(t+\Delta t)\right]$ should be geometrically close to the point at the previous time instant, $\left[w_1(t),\ w_2(t),\ w_3(t)\right]$.

So, the algorithm would look something like this:

Initialize your algorithm by finding $k$, which you said you can do.
Solve for the initial weights $w_1(t), w_2(t), w_3(t)$, which you said you can do.
On subsequent time steps, measure $a(t)$, $b(t)$, and $c(t)$. This defines the plane in $\mathbb{R}^3$ that the filter weight vector can possibly lie upon.
Find the point on the plane that is closest to the filter weights from the previous iteration. Use this point as the new vector of filter weights.
Repeat.

I'm not sure if this will give the desired effect or not (as your inquiry is light on details), but it has an intuitive geometric explanation. It might be worth a try.


Formulas to plot Frequency Impulse Response Coefficients

Assuming you have a set of (windowed) filter coefficients $h[n]$, $n=0,1,\ldots,N-1$, and you want to evaluate the corresponding frequency response at $ M\ge N$ equidistant normalized frequency points

$$f_k=\frac{k}{2M},\quad k=0,1,\ldots,M-1\tag{1}$$

then you need to evaluate the following formula, which is the discrete Fourier transform (DFT):

$$H[k]=\sum_{n=0}^{N-1}h[n]e^{-j\pi nk /M},\quad k=0,1,\ldots, M-1\tag{2}$$

The actual frequencies in Hz are given by $f_k\cdot f_s$, where $f_s$ is the sampling frequency.

Equation (2) is exactly what you get with the Matlab command

[H,w] = freqz(h,1,M)

where $f_k=$ w(k+1)/(2*pi) (because Matlab indices start at $1$, not at $0$).

If you want to evaluate the frequency response at arbitrary, not necessarily equidistant normalized frequencies $f_k$, you need to evaluate

$$H[f_k]=\sum_{n=0}^{N-1}h[n]e^{-j2\pi nf_k},\quad 0\le f_k\le 0.5\tag{3}$$

which for the choice of $f_k$ according to (1) of course reduces to Eq. (2).


Calculating the energy of trapezoid wave

it seems that there power are the same,but there energy are not the same ,because energy=power $\times$ time,and time are different,so their energy are not the same too


It all depends on the transmission medium (the so-called communication channel). Some media do not tolerate components at DC (i.e., zero frequency) or near DC. A typical channel where a carrier frequency greater than zero is required is the radio channel.


What is the relation between kernel functions, kernels used in convolution and null spaces of a matrix?

A kernel in the context of digital signal processing refers to the impulse response $h[n]$ of a filter. Particularly for finite length FIR kernels, filtering is carried out by the convolution operator; hence the kernel of the filter is the kernel of the convolution sum. This kernel maps an input vector $x[n]$ into an output vector $y[n]$ i.e.,
$$ y[n] = h[n] \star x[n] = \sum_{k=-\infty}^{\infty} h[k]x[n-k] ~,$$the sum truncated to a finite length for practical FIR convolutions.

As you have described, a null space of a linear mapping is the set of vectors $x$ whose image is the null vector in the range space of the mapping given by the matrix A. As the name implies it is a set of vectors and not a kernel function (sequence). The only kernel in this context is the matrix A which maps input vectors to output vectors, just as in the context of a convolution operator. $$y = A x$$


CUDA vs FPGA for Image Processing; which is better

Did you have a chance to look at Spiral site?  They try to address this question.  For example, there are special cases where FPGA will run O(N) instead of the customary O(N log N).

Here's a paper that deals with generating FPGA.

All the way at the end of benchmarks page they show performance of a software algorithm enhanced with FPGA-based solvers.

P.S.  It's worth noting that Spiral can be used to optimize algorithms other than FFT, but presented data outlines the popular case only.


How Matched filter is transformed due to frequency modulation?

Normally you wouldn't try to do the match filter directly on the received FM signal.  Typically the receiver path would look something like this-

Received Signal -> Noise/Interference Filter -> FM Demod -> Matched Filter


How to determine stopband of discrete Gaussian, stdev sigma, support N

I found answer at https://courseware.ee.calpoly.edu/~fdepiero/curr/dsp/dsp11.htm . Converted to the present problem, in MATLAB:

kernel2 = fspecial('gaussian', [1 5], 1.3)
spike = [1 0 0 0 0];
[h, w] = freqz(kernel2, spike, 1024);
m = abs(h)
m_db = 10*log10(m/max(m));
plot(w,m_db)

Because the compact support Gaussian had zero crossings, I used the first zero to define the stop band.


Butterworth band-pass filter transfer function

okay, the magnitude frequency response for a general $N$th-order Butterworth lowpass filter is:

$$ \Big|H_\mathrm{LP}(f)\Big|^2 = \frac{1}{1 + \left(\frac{f^2}{W^2} \right)^N} $$

for bandpass filtering, centered at frequency $f_0$, this substitution is made to the LPF:

$$ f \leftarrow \frac{f}{f_0} - \frac{f_0}{f} $$

so then it comes out as:

$$ \Big|H_\mathrm{BP}(f)\Big|^2 = \frac{1}{1 + \left(\frac{\left(\frac{f}{f_0} - \frac{f_0}{f}\right)^2}{W^2} \right)^N} $$

blast that out and you'll have the frequency response for an $2N$th-order Butterworth bandpass filter.


Machine Learning - any suggestions to solve Python rounding errors?

I am not sure what exactly is amiss here (referring to the code) but here is a list of things that might help:

Librosa's mfcc returns a matrix that is n_mfcc x ceil(len(y)/hop_length), so no surprises there.
Divisions between integers can return reals in Python 3 but not Python 2. In any case, if you have a division of 16000/1024 try to enter it as 16000/1024.0 to make sure that what is returned is a real, even if it is going to be rounded in any way later on. In general, work with reals all the way up to the point where you need the rounding and then round at that point.

That's not so important around block_start but it would be around ms_start.

To find the start of a frame containing time offset $to$, you need to $Frn_{start}= \lfloor{\frac{to}{\frac{NFFT-HopLength}{Fs}}}\rfloor$. Where $\lfloor{\cdot}\rfloor$ denotes floor.

Which, yes, it is essentially what you are doing here since int will floor.

Assuming that there is no further mismatch between the hardcoded sampling rates and the actual file that is read from the fisk or a mismatch between the audio file and its text file containing the annotations, I cannot see where is the problem with the labeling.

Can you share a specific example where the frame number is off by 1 or 2?

Hope this helps.

EDIT:

Since the annotation file is a simple text file, can you please collect the files that seem to be having issues and confirm that the annotations are indeed "correct"? To do this you can use Audacity. Audacity accepts simple tab delimited text files with annotation information. The Jamendo corpus includes annotations in space-delimited format and files ending in .lab. The conversion is ultra simple using gnumeric, to import the file and export it as tab delimited. Once you do that, open the MP3 in Audacity and import its labels from the newly exported text file.

In the following screenshot, I have imported the emporte-moi.mp3, converted it to mono and downsampled to 16kHz. This did not seem to have affected the annotations, although certain points, I would have annotated differently because it seems to "cut" when the voice still fades out. But that's a detail probably.

Would be useful to see if there is something going on with the annotations of the files that seem to be "out of sync" with the frame count.


I did it for the DCT-I and DCT-II.
At first I thought it is about circular convolution, but it is not. After desperate attempts to do it by myself I found the article: Convolution Using Discrete Sine and Cosine Transforms (Alternative Link).
In this paper there is a derivation of circular convolution using DTT's, but it is constructed from sum of $C_1^{-1}$ and $S_1^{-1}$, so, it's not what we are interested in.
What I suggest is to look at Fourier transform first, in order to get the idea:
Let we have
$\displaystyle F[t] = \sum_{k=0}^{N-1} f[k]e^{-\frac{i2π}{N}tk}$,
$\displaystyle G[t] = \sum_{k=0}^{N-1} g[k]e^{-\frac{i2π}{N}tk}$.
So then
$$\displaystyle ℱ\{f*g\}[t] = \sum_{x=0}^{N-1} f[x]e^{-\frac{i2π}{N}tx}\cdot\sum_{y=0}^{N-1} g[y]e^{-\frac{i2π}{N}ty} = \sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f[x]g[y]e^{-\frac{i2π}{N}t(x+y)}$$
Now, let $x+y = m$, $m = 0\ldots2N-2$ and then group all $(x,y)$ which share the same $m$. So you will get:
$$\displaystyle \sum_{m=0}^{2N-2}\biggl(\sum_{0\leq x,y\leq N-1\\x+y=m}f[x]g[y] \biggl)   e^{-\frac{i2π}{N}tm} \tag{1}$$
Because $x, y$ are both positive and less than $N$, then sum inside is just linear convolution at $m$. Let me denote it as $A[m]$, then everything is going to be rewritten in that form
$$\sum_{m=0}^{2N-2} A[m] e^{-\frac{i2π}{N}tm} \tag{2}$$
DFT kernel is N-periodic that is why we can group terms another time:
$$\displaystyle \sum_{m=0}^{N-1} (A[m] + A[m+N])e^{-\frac{i2π}{N}tm}$$
$A[m]$ captures $A[0]\ldots A[N-1]$, $A[m+N]$ captures $A[N]\ldots A[2N-1]$. But remember that $A[m]$ is linear convolution of two sequences of size $N$ and therefore size is $N + N - 1$ then $A[2N-1] = 0$. So these two sums are equivalent. Also it is a DFT of $A[m] + A[m+N]$. What is it? Right, circular convolution, as expected.
Now plug in DHT kernel instead of DFT into $(1)$ and you will get the same result ─ circular convolution. That is because DHT and DFT kernels are quite similar.
Well, if you have 2 signals $f$, $g$ with lengths $Q$ and $P$, and $A = f*g$, with $Q+P-1 = L$. Now, add zeros to $f$ up to $L$, do the same for $g$. Now $f$, $g$, and $A$ are $L$ periodic signals, so we can take the DFT of them. DFT of $A$ is going to be:
$$\sum_{m=0}^{L-1} A[m] e^{-\frac{i2π}{L}tm} = \sum_{m=0}^{L-1}\biggl(\sum_{0\leq x,y\leq L-1\\x+y=m}f[x]g[y] \biggl)   e^{-\frac{i2π}{L}tm} = \sum_{x=0}^{L-1}\sum_{y=0}^{L-1}f[x]g[y]e^{-\frac{i2π}{L}t(x+y)} \tag{3}$$
So, via $(2)$ we can get the rule and expression for the convolution and $(3)$ shows that zero-padding up to convolution length is required to get linear convolution using the same rule. We're going to apply it to the DCT's and DST's.
Because there are many varieties of DCT's and DST's, I will mention what I used:
DCT-I (N+1) points:
$\displaystyle C_1\{f\}[n] = \sum_{k=0}^{N} c_k f[k]\\cos\Bigl[\frac{π}{N}nk\Bigl], \text{where}\
c_k =
\left\{
\begin{array}{l}
\tfrac12, \qquad &k (\mod N)=0, \\
1,\qquad &\text{otherwise} \\
\end{array}
\right. \tag{4}$
DST-I (N+1) points:
$$S_1\{f\}[n] = \sum_{k=0}^{N} f[k]\sin\Bigl[\frac{π}{N}nk\Bigl]=\sum_{k=0}^{N} c_k f[k] \sin\Bigl[\frac{π}{N}nk\Bigl] \tag{5}$$
(N-1) for data: $f[1] \ldots f[N-1]$ it doesn't matter what is placed in $f[0]$ and $f[N]$.
DCT-II (N+1) points:
$$C_2\{f\}[n] = \sum_{k=0}^{N} f[k]\cos\Bigl[\frac{π}{N}\Bigl(k+\tfrac12\Bigl)n\Bigl] \tag{6}$$
N for data: $f[0]\ldots f[N-1]$ and $f[N] = 0$
DST-II (N+1) points, N for data:
$f[0]\ldots f[N-1]$ and $f[N] = 0$:
$$S_2\{f\}[n] = \sum_{k=0}^{N} f[k]\sin\Bigl[\frac{π}{N}\Bigl(k+\tfrac12\Bigl)n\Bigl] \tag{7}$$
Inverse DCT-I:
$$C_1^{-1} = \frac{2}{N}C_1 \tag{8}$$
Inverse DST-I:
$$S_1^{-1} = \frac{2}{N}S_1 \tag{9}$$
Inverse DCT-II (N+1) points:
$$C_2^{-1}\{F_{C_2}\}[k] = \frac{2}{N}\sum_{n=0}^{N} c_n F_{C_2}[n]\cos\Bigl[\frac{π}{N}n\Bigl(k+\tfrac12\Bigl)\Bigl] \tag{10}$$
$C_2^{-1}\{C_2\{f\}\}[k] = f[k]$ except for $k=N$ there it is going to be $f[N-1]$.
Inverse DST-II (N+1) points:
$$C_2^{-1}\{F_{S_2}\}[k] = \frac{2}{N}\sum_{n=0}^{N} c_n F_{S_2}[n]\sin\Bigl[\frac{π}{N}n\Bigl(k+\tfrac12\Bigl)\Bigl] \tag{11}$$
$S_2^{-1}\{S_2\{f\}\}[k] = f[k]$ except for $k=N$ there it is going to be $-f[N-1]$. Also, they are not interchangeable.
Everything mentioned above is just Wikipedia ones, but modified to achieve more similar look for easier coding and derivation.
So, now again, applying the rule $(2)$ for the DCT-I we get:
$$\sum_{m=0}^{2N} A[m] c_m \cos\Bigl[\frac{π}{N}nm\Bigl] = \sum_{m=0}^{N} p_m(A[m] + A[2N-m]) c_m \cos\Bigl[\frac{π}{N}nm\Bigl]$$
$p_m = \tfrac12$ if $m = N$ else $p_m = 1$
from the other hand we have
$$\sum_{m=0}^{2N} \frac{A[m]}{c_m} c_m \cos\Bigl[\frac{π}{N}nm\Bigl] = \sum_{x=0}^{N}\sum_{y=0}^{N} f[x]g[y]\cos\Bigl[\frac{π}{N}n(x+y)\Bigl]$$
$$\sum_{x=0}^{N}\sum_{y=0}^{N} f[x]g[y]\cos\Bigl[\frac{π}{N}n(x+y)\Bigl] = \\\sum_{x=0}^{N}\sum_{y=0}^{N} f[x]g[y]\biggl(\cos\Bigl[\frac{π}{N}nx\Bigl]\cos\Bigl[\frac{π}{N}ny\Bigl] - \sin\Bigl[\frac{π}{N}nx\Bigl]\sin\Bigl[\frac{π}{N}ny\Bigl]\biggl) = \\ \biggl(\sum_{x=0}^{N} f[x]\cos\Bigl[\frac{π}{N}nx\Bigl]\biggl)\cdot\biggl(\sum_{y=0}^{N} g[y]\cos\Bigl[\frac{π}{N}ny\Bigl]\biggl) - \biggl(\sum_{x=0}^{N} f[x]\sin \Bigl[\frac{π}{N}nx\Bigl]\biggl)\cdot\biggl(\sum_{y=0}^{N} g[y]\sin\Bigl[\frac{π}{N}ny\Bigl]\biggl)$$
It is now easy to see that it can be rewriten in a such way:
$$C_1\{\frac{f[m]}{c_m}\}[n]\cdot C_1\{\frac{g[m]}{c_m}\}[n] - S_1\{f[m]\}[n]\cdot S_1\{g[m]\}[n]$$
Let $\frac{f[m]}{c_m} = f^{/}$ then we can write everything in this form:
$$\frac{p_m}{c_m}(A[m] + A[2N-m]) ⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}$$
alternative version
$$δ_{0m}(A[0] + A[2N]) + (A[m] + A[2N-m]) ⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}$$
or even
$$(A[m] + A[2N-m]) ⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1} - \frac{f[0]g[0] + f[N]g[N]}{2}$$
Now, remember the $(3)$? If you add zeroes, so the size of f and g is equal to the convolution's length, then you have:
$$(f*g)^{/}⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}$$
It looks a little bit different to what is presented in the paper, but it is because I used different version of DCT.
For DCT-II we have the rule (remember that $f[N] = g[N] = 0$):
$$\sum_{m=0}^{2N} A[m] \cos\Bigl[\frac{π}{N}\Bigl(m+\tfrac12\Bigl)n\Bigl] = \sum_{m=0}^{2N-2} A[m] \cos\Bigl[\frac{π}{N}\Bigl(m+\tfrac12\Bigl)n\Bigl]  = \sum_{m=0}^{N-1} (A[m] + A[2N-1-m]) \cos\Bigl[\frac{π}{N}\Bigl(m+\tfrac12\Bigl)n\Bigl] $$
And the expression is going to be:
$$\sum_{m=0}^{2N} A[m] \cos\Bigl[\frac{π}{N}\Bigl(m+\tfrac12\Bigl)n\Bigl] =
\sum_{m=0}^{2N} A[m] \Bigl(\cos\Bigl[\frac{πn}{2N}\Bigl]\cos\Bigl[\frac{π}{N}kn\Bigl] - \sin\Bigl[\frac{πn}{2N}\Bigl]\sin\Bigl[\frac{π}{N}kn\Bigl]\Bigl) = \\
\cos\Bigl[\frac{πn}{2N}\Bigl](F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}) - \sin\Bigl[\frac{πn}{2N}\Bigl](F^{/}_{C_1}G_{S_1} + F_{S_1}G^{/}_{C_1}) = \\
=F^{/}_{C_1}(\cos\Bigl[\frac{πn}{2N}\Bigl]G^{/}_{C_1} - \sin\Bigl[\frac{πn}{2N}\Bigl]G_{S_1}) - F_{S_1}(\cos\Bigl[\frac{πn}{2N}\Bigl]G_{S_1}+ \sin\Bigl[\frac{πn}{2N}\Bigl]G^{/}_{C_1}) = \\
=F^{/}_{C_1}G_{C_2} - F_{S_1}G_{S_2} = G^{/}_{C_1}F_{C_2} - G_{S_1}F_{S_2} $$
Now we finally can write it as:
$$A[m] + A[2N-1-m] ⟷ F^{/}_{C_1}G_{C_2} - F_{S_1}G_{S_2}$$
And again, if you add enough zeroes, then it will be:
$$f*g ⟷ F^{/}_{C_1}G_{C_2} - F_{S_1}G_{S_2}$$
You can play with it here: https://repl.it/@dogtyping/dct#main.py — to make sure that it does work.
Thanks for the question, it was a lot of fun.


As you've indicated there may be some mismatching in the filter bands resulting in undesirable results. If one of the filters is heavily attenuating the frequencies of the other filter, there's little you can do to reconcile this by just cascading the filters together.

Given you'd like to maintain the response of both filters independent of the other, you need to add the filter outputs together.

Here's some code that shows the result of one filter heavily attenuating the frequencies of the other. It also shows the equivalence of cascading two filters either by connecting the output to the input of the next stage or just convolving the filter taps.

b1 = fir2(2048, [0 0.1 0.15 0.3 0.35 1], [0 0 1 1 0 0]);
b2 = fir2(128, [0 0.4 0.45 0.6 0.65 1], [0 0 1 1 0 0]);

x = randn(1, 32768);

% Parallel sum
y1 = filter(b1, 1, x) + filter(b2, 1, x);

% Cascade
y2 = filter(b2, 1, filter(b1, 1, x));

% Equivalent cascade
y3 = filter(conv(b1, b2), 1, x);

figure, hold all
subplot(311), plot(20*log10(abs(fft(y1))))
subplot(312), plot(20*log10(abs(fft(y2))))
subplot(313), plot(20*log10(abs(fft(y3))))


There is nothing wrong with your calculation. However, by definition, $x[n]=0$ for $n<0$, because otherwise we couldn't have used the unilateral $\mathcal{Z}$-transform in the first place. So you get

$$\mathcal{Z}\big\{x[n-k]\big\}=x^{-k}X(z),\quad k> 0\tag{1}$$

Eq. $(1)$ is valid for the unilateral as well as for the bilateral $\mathcal{Z}$-transform.

Things are different for a time advance. The unilateral $\mathcal{Z}$-transform of $x[n+k]$, $k>0$, is

$$\begin{align}\mathcal{Z}\big\{x[n+k]\big\}&=\sum_{n=0}^{\infty}x[n+k]z^{-n}\\&=\sum_{n=k}^{\infty}x[n]z^{-(n-k)}\\&=z^k\left(X(z)-\sum_{n=0}^{k-1}x[n]z^{-n}\right),\quad k>0\tag{2}\end{align}$$

For the bilateral $\mathcal{Z}$-transform we simply have

$$\mathcal{Z}\big\{x[n+k]\big\}=z^kX(z)\tag{3}$$


How to implement IDFT using the positive frequencies of DFT, modulation and resampling on MATLAB?

Here is the IDFT with normal and positive only frequencies if it helps:

N=4;
A=rand(N-1,1) + j*rand(N-1,1);   % Random positive carrier
X=[0; A; 0; flipud(conj(A))];    % Hermitian symmetric DFT vector
Xp=[0; 2*A; 0; zeros(size(A))];  % Positive freq only DFT vector

x = ifft(X);         % Normal IDFT
xa = ifft(Xp);       % Analytic signal
xp = real(ifft(Xp)); % Same as normal IDFT


How do I calculate the eigenvalues of a covariance matrix which contains harmonic functions

The included document does not say $ce^{-\frac{\phi_i-\phi_j}{\rho}}$ has the form of $e^{-in\phi_i}$. The included document uses $C$ for three completely different things:

$C_{ij}$ is the correlations between $i$ and $j$ items.
$C(\phi_i-\phi_j)$ a component of that correlation, dependent on the angles between the two items $i$ and $j$.
$C_n$ the eigenvalues of the covariance matrix made up of the $C_{ij}$.

The term $e^{-in\phi_i}$ is a component of $C_n$ the eigenvalues / Fourier components expression in (31).

The idea would be to form the matrix $C_{ij}$ and then take the eigensystem decomposition of that matrix.


How do you interpret FFT of an FFT of a discrete signal?

Essentially, what you do is performing the DFT twice. Performing DFT twice amounts to time-reversing the signal:

N = 16; sN= np.sqrt(N)
x = np.random.randn(N)

X = np.fft.fft(x) / sN

x2 = np.fft.fft(X) / sN

plt.plot(x.real)
plt.plot(x2.real)

blue is original signal, green is two-times DFT. As you can see, the signals are time-reversed in the sense of $y[n]=x[-n \mod N]=x[N-n], n=0...N-1$. (y is green, x is blue). Here, I used the unitary version of the DFT, to keep the signals in the same scale. Otherwise, with the "normal" definition of DFT, the green signal would be scaled by 16.


Does Discrete Wavelet Consider Sampling Frequency or Cutoff Frequency of Bandpass Filter?

THe DWT does not really know about the actual sampling. It do cares about the relative frequency span. And it has not knowledge about pre-processing, like linear filtering.
So in your case, the detail coefficients of the first level of wavelet coefficients correspond,  more or less, to the [128/2-128]/2 frequency range.
Because of the initial [0-30] low-pass filter, the above [32-64] first wavelet subband range is likely to contain little useful information.


Since the data records $x_i$ are uncorrelated realizations of the same random process, $\hat P^{per}_{(i)}(e^{j\omega})$ are all uncorrelated random variables with identical means and variances (given by the mean and variance of the Bartlett Method's PSD estimator). So, for any fixed $\omega$, $\mathbf E[\hat P^{per}_{(i)}(e^{j\omega})]$ are equal $\forall i$ and also $\mathrm{Var}[\hat P^{per}_{(i)}(e^{j\omega})]$ are equal $\forall i$. So,
\begin{eqnarray}
\mathbf E[\hat P_{x}(e^{j\omega})] &=& \mathbf E\left[\frac{1}{K} \sum_{i=1}^K\hat P^{per}_{(i)}(e^{j\omega})\right] \\
&=& \frac{1}{K}\sum_{i=1}^K \mathbf E\left[\hat P^{per}_{(i)}(e^{j\omega})\right] \\
&=& \frac{1}{K} K \;\mathbf E\left[\hat P^{per}_{(i)}(e^{j\omega})\right]\\
&=& \mathbf E\left[\hat P^{per}_{(i)}(e^{j\omega})\right]
\end{eqnarray}

and similarly,

\begin{eqnarray}
\mathrm{Var}[\hat P_{x}(e^{j\omega})] &=& \mathrm{Var}\left[\frac{1}{K} \sum_{i=1}^K\hat P^{per}_{(i)}(e^{j\omega})\right] \\
&=& \frac{1}{K^2}\sum_{i=1}^K \mathrm{Var}\left[\hat P^{per}_{(i)}(e^{j\omega})\right] \\
&=& \frac{1}{K^2} K \;\mathrm{Var}\left[\hat P^{per}_{(i)}(e^{j\omega})\right]\\
&=& \frac{1}{K}\mathrm{Var}\left[\hat P^{per}_{(i)}(e^{j\omega})\right].
\end{eqnarray}


It's simpler than you think. When we discretize frequencies, we get frequency bins. So, when you discretize your Fourier Transform: $$e^{-j\omega} \rightarrow e^{-j{2\pi k}/{N}}$$ Our continuous frequencies become $N$ discrete bins.

This is exactly why the following is true: $$n^{th}bin = n*\dfrac{sampleFreq}{num(DFT points)}$$ Note that the FFT represents frequencies 0 to sampleFreq Hz.


Performed the Danielson-Lanczos shuffling for FFT, but I don't know what to do next

I am not so sure about what you're asking for, but I will try to tell you something about this point of your question:

  The problem is that I'm a little confused on what to do with the resultant vector.

You have not a resultant vector, you have two vectors [0,2,4,6], [1,3,5,7] and by the D-L lemma you can compute the  DFT of the original vector as a linear combination of the DFTs of these two vectors.

Taking into account that the two shorter DFTs are periodic and that the combination coefficients for the second half are the negative of the coefficients used for the first half, we are lead either to a very simple recursive implementation of the FFT algorithm or to more complicated implementations that don't use recursion. Let's stay simple...

In a pythonesque pseudocode, the recursive implementation can be written as simply as

def fft(x,n,e):
    if n==1:
        X = x # the DFT of a scalar is the scalar itself
        return X
    else:
        X0, X1 = fft(x[0::2],n/2,e[0::2]), fft(x[1::2],n/2,e[0::2])
        X0, X1 = X0 + e*X1, X0-e*X1
        X = concatenate(X0,X1)
        return X

where

x is the vector of which we want to compute the DFT
n is its length
e is a vector with the so called _twiddle_factors_, that is one half of the n-th roots of unity (from the lower part of the unit circle for forward DFT, from the upper otherwise)

and the [start:stop:step] notation stands for a selection of the elements of a vector, with a missing stop meaning to the end.

To call fft you have to provide the initial vector of twiddle factors

X = fft( x,  n, exp(-2*pi*1j*(0..n/2-1)/n)) # forward
x = fft(X/n, n, exp(+2*pi*1j*(0..n/2-1)/n)) # inverse

where 1j id the imaginary unit and (0..n/2-1) indicates an aritmetic progression from 0 to n/2-1, extremes included.


Question on covariance matrix of 2 spatial signals

If you have two signal vectors $x_1[n]$ and $x_2[n]$ each of $N$ elements, then there are two different things we can consider.

How do the quantities $\displaystyle\sum_{n=1}^N x_i[n]x_j[n],~ i, j \in \{1,2\}$ compare?  In particular, when the signals are noisy and the noises can be considered
to be jointly stationary (or jointly wide-sense stationary), these quantities can be used to estimate the noise variances in the two signals as well as the covariance of the noises at any fixed sampling time.  This is what you get from the $2\times 2$
covariance matrix
$$R_{2\times 2} = \left[\begin{matrix} \sigma_1^2 & C\\
C & \sigma_2^2\end{matrix}\right].$$
The noise in $x_1[n]$ has variance
$\sigma_1^2 = R_{1,1}$ which might be different
from $R_{2,2} = \sigma_2^2$, the variance of the noise in $x_2[n]$.  But the noises are
correlated with covariance $R_{1.2}=R_{2,1} = C$.  Now, if we plan on doing things
with just what happens at $n$, ignoring whatever might be happening at $n-1$ or $n+1$
etc., then this is all the information we need.
Unless the noise is known to be (or assumed to be) white noise so that noise
samples from different sampling instants are independent (and therefore uncorrelated)
or we simply assume uncorrelated noise samples, there is information that we
are ignoring by not considering the correlation between $x_1[n]$ and $x_1[m]$,
samples from the same process at different times or locations, and the correlation between $x_1[n]$ and $x_2[m]$, samples from the two processes at different times
or locations.  This additional information might lead to a better estimate/solution.
Now  we have a total of $2N$ noise samples,
and therefore a $2N \times 2N$ covariance matrix to consider.  If we
arrange matters the way the authors did, we have
$R_{\text{full}} = E[XX^T]$ where
$$X = (x_1[1],x_1[2],\ldots,x_1[N],x_2[1],x_2[2],\ldots,x_2[N])^T
= (\mathbf x_1,\mathbf x_2)^T$$
and so
$$R_{\text{full}} = \left[
\begin{matrix}R_{\mathbf x_1,\mathbf x_1} & R_{\mathbf x_1,\mathbf x_2}\\
R_{\mathbf x_2,\mathbf x_1} & R_{\mathbf x_2,\mathbf x_2}\end{matrix}
\right]$$ where $R_{\mathbf x_i,\mathbf x_j} = E[\mathbf x_i\mathbf x_j^T]$.
Note that $R_{\mathbf x_i,\mathbf x_j}$ is, in essence, the
cross-correlation function of $(x_i[1],x_i[2],\ldots,x_i[N])$
and $(x_j[1],x_j[2],\ldots,x_j[N])$ if $i \neq j$ and the
autocorrelation function if $i = j$.  If the noise processes are white and
uncorrelated except when $n = m$, then
$$R_{\text{full}} \rightarrow R_{\text{simple}} = \left[ \begin{matrix}\sigma_1^2I & CI\\
CI & \sigma_2^2I \end{matrix}\right]$$
where $I$ is the $N\times N$ identity matrix, and
$\sigma_1^2, \sigma_2^2$ and $C$ are as defined in Item 1 above.
How realistic this noise model might be is something for the end user
to determine. If the model is realistic, then nothing is gained
by looking at the $2N\times 2N$ matrix $R_{\text{full}}$
since all the information is there
in the $2\times 2$ matrix $R_{2\times 2}$ of Item 1 above.  Ditto if the model is
unrealistic but we don't intend to (or are unable to)
use all the information in the full
$2N\times 2N$ matrix $R_{\text{full}}$; we will make do
with just $\sigma_1^2, \sigma_2^2$ and $C$ of Part 1
for which we don't need $R_{\text{full}}$ or $R_{\text{simple}}$,
just $R_{2\times 2}$.


Use Scale Space Representation to Filter Single Image

I used factor 5-6 going from the Standard Deviation (Std) of the Gaussian to the Kernel:
radius = ceil(6 * STD);

Though it means more computation power is required.


It looks like you are mixing up two techniques: MinBLEP and band-limited wavetables.

Band-limited wavetables

Let's assume you have a 256 sample long array with a waveform of interest - for example a naive sawtooth, with the first sample in your array at -1 and the last sample at 1.

To replay this data at a variable pitch, two techniques are possible:

First, you can control the rate at which the data is played. For example, to get a 440 Hz pitch, you replay the 256 samples verbatim, at a rate of 112640 Hz. This is possible with custom hardware, and is the technique used in early wavetable synthesizers (PPG Wave, Prophet VS). This requires one clock/DAC per oscillator, and is rather expensive!

Or you use a fixed sample-rate, and you modify the pitch by interpolation through the wavetable. Let's say the sample rate is 48kHz. Playing the data one sample after the other yields a tone with a frequency of $48000 / 256 = 187.5$ Hz. To play a tone at 440Hz, you would need to play the data at $440 / (48000 / 256)=2.35$ its original rate, and that's where you get into trouble. Two ways of understanding why:

You can consider the discontinuity. You'll reach the end of the table and "wrap around" between the 109th and 110th sample. So the discontinuity should actually be located between these samples.
Or you can think in terms of frequency content: by playing the data at 2.35 its original rate, any harmonic it contains at a frequency above nyquist / 2.35 will be aliased.

There are two ways of dealing with the problem. A first approach is to use band-limited interpolation, which is a way of replaying the data in such a way that aliased components are filtered. Another solution is to pre-compute several variants of your waveform in which the highest harmonics are progressively filtered. For example, when you have to playback the data at twice its original rate, you use a different sawtooth wavetable in which all harmonics above nyquist/4 are filtered-out. Obviously, you can't have a table ready for every possible playback rate, so a common solution is to pre-compute a dozen of such tables (for example one per octave), and interpolate between them. Korg did that on their DW series, I used it on my DIY AVR-based synths.

There are different ways of building such a band-limited table. You can work by additive synthesis (adding the right number of harmonics), or try to derive the band-limited waveforms by operations on a band-limited impulse (sinc waveform).

minBlep

This technique is completely unrelated to wavetable synthesis. The idea is to precompute a table containing a band-limited and causal step waveform. This computation is done once and for all, and involves an FFT/IFFT - that's probably where you encountered these operations, which are used only at this offline step. The code you link to in your question only covers this step!

Then, in real-time, you generate the naive waveform, and whenever there is a discontinuity in the waveform, you add, over the next 40 samples or so, samples from the band-limited step table to smooth the edge. Which samples you need to add from the minBlep table are determined by the fractional sample at which the discontinuity occurs. No look-ahead is involved: the point of minBLEP is that the correction is causal, so you don't have to start adding it before the discontinuity occurs.

Note that the computational cost of this method will increase as the frequency goes up, since the discontinuities will be more frequent - you'll have to keep track of an increasing number of overlapping corrections.

I wouldn't recommend minBLEP by the way. PolyBLEP uses a polynomial approximation of the band-limited step function and yields very good results. No lookup data is required and the code is easier to write because the correction only lasts for a handful of samples.

The two together?

It would be a strange idea to mix the two techniques together, but someone could use minBLEP offline to precompute a band-limited sawtooth for various band-limiting ratios. For example, a 256 sample-long naive sawtooth, then a 256 sample-long sawtooth with no harmonic above 0.25, then a 256 sample-long sawtooth with no harmonic above 0.125... And then one could use wavetable playback to replay them - selecting/crossfading one of these tables depending on the required playback rate. This would be a bit odd, there are simpler ways of building wavetables!


Phase Correlation - Poor Performance on Noisy/Blurred Images?

One-dimensional version

The one-dimensional version that you list won't work. When there is a large enough shift in images (more than one or two pixels in real-world images), there will be nothing relating the column pixels.

For an example of this, try:

I5 = rand(100,100)*255;
I6 = zeros(100,100);
I6(11:100,22:100) = I5(1:90,1:79);

So that we have I5:

and I6:

Then the one-dimensional phase correlation is just:

whereas the two-dimensional phase correlation is:

It's a bit hard to see, but there is a very high peak in the bottom-righthand corner of the image. No clear peak exists in the one dimensional version.

Why doesn't smoothing help? #1

What correlation is trying to do is to find "similar" variations in each image.  If the underlying signals are sufficiently random, then this will work well: the correlation of white noise with itself gives a really nice peak at the origin, and close-to-zero elsewhere.

Smoothing a "random" image with a Gaussian will have the effect of smoothing out the correlation your expecting --- spreading the energy in any peaks over a wider area.

Smoothing has the opposite effect of "pre-whitening" the image. Pre-whitening (as the name suggests) tries to make the image more like white noise --- which has the best form if we are doing correlation-based detection (in that the peak is well-localized).

What you are better off doing is using the matlab diff operation for a simplistic, but surprisingly effective, way of pre-whitening images.

See this example.

Why doesn't smoothing help? #2

Why does smoothing cause the extra peaks?

If you smooth each image with a kernel $k(x,y)$ then we get:
$$
h_a = g_a \star k\\
h_b = g_b \star k
$$
where $\star$ is convolution.

Now,
$$
H_a = K G_a\\
H_b = K G_b\\
R = \frac{H_a H_b^*}{|H_a H_b^*|} = \frac{|K|^2G_a G_b^*}{|K|^2|G_a G_b^*|}\\
=  \frac{G_a G_b^*}{|G_a G_b^*|}
$$

What I suspect is happening (though I am not sure) is that perhaps your kernel has near-zero values in the frequency domain, causing numerical problems?

If I apply a kernel:

K = one(5,5);

to my random images, then I get:

for the two-dimensional correlation, which makes the peak more spread out, but does not exhibit the problems you are seeing.


Spectral efficiency of with OFDM vs without OFDM?

The spectral efficiency of OFDM is strictly worse than that of a pulse-shaped QAM signal with the same rate.

OFDM requires a guard interval, on which no useful information is transmitted.

It is common to dedicate several subcarriers to pilot signals.

It is often not pulse-shaped (or, more accurately, pulse-shaped with a rectangular pulse).

Since each non-pilot subcarrier is essentially a narrowband QAM signal, you can see that its spectral efficiency cannot be better. However, the point of OFDM is to allow for low-complexity frequency-domain equalization, with spectral efficiency given a less important role.


It looks to me as if you are asking two completely separate questions here.

First question: Indeed the two formulas aren't identical; it isn't obvious that the Goertzel algorithm computes what it does. You can find a sketchy proof based on the $z$-transform on the Wikipedia page about the Goertzel algorithm, or you can do a more direct proof by mathematical induction. Or you can just take it on trust -- for most practical purposes all that matters is that it does what it does, rather than why.

Second question:

With both the Goertzel algorithm and the complete DFT, |X|^2 is the power at a particular frequency. This is a sensible thing to compute because, e.g., when you do a DFT the sum of all the powers (|Fourier coefficients|^2) equals the total power in the signal (= sum of |signal values|^2, up to a normalizing constant factor).

If the frequency you're using is one of the ones processed by a DFT, the output of the Goertzel algorithm is (aside from "numerical" considerations like roundoff error) exactly the same as one bin of the DFT.


Understanding the Parameters of the Bilateral Filter

A small discussion of the details you've asked for:

The bilateral ﬁlter is controlled by two parameters: $\sigma$d and $\sigma$r.

• As the range parameter $\sigma$r increases, the bilateral ﬁlter becomes closer to Gaussian blur because the range Gaussian is ﬂatter i.e., almost a constant over the intensity interval covered by the image.

• Increasing the spatial parameter $\sigma$d smooths larger features.

An important characteristic of bilateral ﬁltering is that the weights are multiplied, which
implies that as soon as one of the weight is close to 0, no smoothing occurs. As an example, a large spatial Gaussian coupled with narrow range Gaussian achieves a limited smoothing although the ﬁlter has large spatial extent. The range weight enforces a strict preservation of the contours.

"The Norm or the double bars indicate the gaussian distance in the equation.This distance is deﬁned by Gσ(||p − q||), where σ is a parameter deﬁning the extension of the neighborhood."

Since (i,j) and (k,l) are simply spatial points on an image, they will vary from the start of the image to the end of the image. Typically this would be dependent on your indexing measure, for example in a matrix, you've have to range from 0 to 5 for a 5x5 image. The actual values of k,l matter little. Its the relative value which matters far more.


this is John BG

1.- You have to work with Signal Power, not just input or output:

Circuit resonance is measured with the quality factor Q and does not only have to do with the signal you send in, or the signal that bounces back, but with those AND the impedance that such signals find.

Literature reference Microwave Engineering (Pozar) 1 defines Circuit Resonance the following way

2.- Measure the impedance:

Without the impedance, there's no way you can assess Q.

You are focusing your effort on the test signal and it's spectrum, but there's not a single reference to the impedance of the circuit you are measuring.

There's no way you can measure Q without at least an additional signal, than can be the impedance that the test input faces, or alternatively the reflected signal, that again you haven't mentioned.

3.- Narrow band or wide band? doesn't really matter

Despite the definition of resonance is done with a central frequency f0, it may be the case that the circuit resonates on each of the frequencies within the Chirp signal you use to poke the system.

4.- Use Chirp with constant amplitude

From the amplitude of the signal you shown in the question, the amplitude decreases with increasing frequency, that complicates calculations.

5.- Wrapping Suggestion:

Why don't you

start with a constant amplitude input chirp,
measure the return signal
from there get the impedance
calculate the 'sharpness' of the impedance notch over frequency,
which is actually what defines how resonant is a circuit,
from there you can tell whether the circuit under test is resonant
throughout the chirp signal frequency band or only on certain
frequencies.
Once you have all this working, you may want to start playing with
the amplitude of the chirp input signal.

If you follow the above steps, please let us know the results you obtain updating this QA with further comments, whether it worked for you, or if not what code have you written on such regard.

If you send me by email the MATLAB code you used to generate the question signal, I will have a look and look for a way to get the reflected signal.

John BG

Literature Reference (1):

https://www.amazon.co.uk/Microwave-Engineering-Written-Publisher-Hardcover/dp/B00SLRKBV2/ref=sr_1_1?ie=UTF8&qid=1529281040&sr=8-1&keywords=David+Pozar+Microwave+Engineering

Excerpts of the book can be read here: http://www2.electron.frba.utn.edu.ar/~jcecconi/Bibliografia/Ocultos/Libros/Microwave_Engineering_David_M_Pozar_4ed_Wiley_2012.pdf

There's also a solutions manual available here: https://www.scribd.com/doc/176505749/Microwave-engineering-pozar-4th-Ed-solutions-manual


Understanding noise removal method using wavelets

Your signal (with initial par x0 =0.1) is already noise like and high frequency. It will be hard to distinguish it from the added white noise... One thing you can do is to interpolate (resample) the time series by a large enough factor and then later add the white noise. This will artifically help to separate the noise spectrum and your signal spectrum but the signal lengths will also be increased. whether it is what you have to do is up to you !
The following modification apparently improves the noise removal, but fundamentlaly it's separating the noise spectrum from the signal. So whether this is a viable option is upto your applications:
 M = 50;        % number of data points (length of the time series)

 x(1) = 0.5;    % initial condition (can be anything from 0 to 1)
 for n = 2:M,   % iterate
     x(n) = 4*x(n-1)*(1-x(n-1));
 end

 U = 10;       % interpolation factor
 xU = resample(x,U,1);  % just interpolate the obtained sequence

 % add noise onto the interpolated sequence xU
 x_noise = awgn(xU , 10 , 'measured');

 %denoise using wavelet
 denoised = wdenoise(x_noise, 3,'Wavelet','db3',...
     'DenoisingMethod','Bayes',...
     'ThresholdRule','Median',...
     'NoiseEstimate','LevelIndependent');

 denoised = resample(denoised,1,U);   % downsample de-noised sequence back

 figure
 plot(x_noise(1:10:end))   % down-sample noisy seqeunce on the fly for displaying
 axis tight
 hold on
 plot(denoised,'r')
 plot(x,'c--');
 legend('noisy','denoised','clean');

 fb = cwtfilterbank('SignalLength',M);
 [cfs,frq] = wt(fb,denoised);

The result looks like :


What Is the Point of Doing the Zero Padding?

You may think of it as efficient way to apply Dirichlet Window Based interpolation in the Fourier Domain.

The advantage of applying the interpolation using Zero Padding in the Time Domain is very simple - it is simpler and more computationally efficient.


As the returned value shows, the first spectrum is computed at time $t=0.15$ s. With a sampling rate of 1000 Hz, that means that the first analysis window is centred at sample $n=150$. So the first analysis frame spans the samples $n \in [0,300[$.

Because your overlap is set to 25 samples, the second analysis frame spans the samples $n \in [275,575[$, thus overlapping in its first 25 samples with the first analysis frame. The second analysis frame is then centred at time $t=0.425$ s or in samples $n=425$.

The third analysis frame again overlaps 25 samples with the second one. Thus, it spans the samples $n \in [550,850[$ and is centred at time $t=0.7$ s or in samples $n=700$.

The next analysis frame would span the samples $n \in [825,1125]$ and is not computed because it would go over the range of samples in the input signal.


One important thing to understand is that after extracting the keypoints, you only obtain information about their position, and sometimes their coverage area (usually approximated by a circle or ellipse) in the image. While the information about keypoint position might sometimes be useful, it does not say much about the keypoints themselves.

Depending on the algorithm used to extract keypoint (SIFT, Harris corners, MSER), you will know some general characteristics of the extracted keypoints (e.g. they are centered around blobs, edges, prominent corners...) but you will not know how different or similar one keypoint is to the other.

Here's two simple examples where only the position and keypoint area will not help us:

If you have an image A (of a bear on a white background), and another image B, exact copy of A but translated for a few pixels: the extracted keypoints will be the same (on the same part of that bear). Those two images should be recognized as same, or similar.

But, if the only information we have is their position, and that changed because of the translation, you can not compare the images.
If you have an image A (let's say, of a duck this time), and another image B, exactly the same duck as in A except twice the size: the extracted keypoints will be the same (same parts of the duck). Those are also same (similar) images.

But all their sizes (areas) will be different: all the keypoints from the image B will be twice the size of those from image A.

So, here come descriptors: they are the way to compare the keypoints. They summarize, in vector format (of constant length) some characteristics about the keypoints. For example, it could be their intensity in the direction of their most pronounced orientation. It's assigning a numerical description to the area of the image the keypoint refers to.

Some important things for descriptors are:

they should be independent of keypoint position

If the same keypoint is extracted at different positions (e.g. because of translation) the descriptor should be the same.
they should be robust against image transformations

Some examples are changes of contrast (e.g. image of the same place during a sunny and cloudy day) and changes of perspective (image of a building from center-right and center-left, we would still like to recognize it as a same building).

Of course, no descriptor is completely robust against all transformations (nor against any single one if it is strong, e.g. big change in perspective).

Different descriptors are designed to be robust against different transformations which is sometimes opposed to the speed it takes to calculate them.
they should be scale independent

The descriptors should take scale in to account. If the "prominent" part of the one keypoint is a vertical line of 10px (inside a circular area with radius of 8px), and the prominent part of another a vertical line of 5px (inside a circular area with radius of 4px) -- these keypoints should be assigned similar descriptors.

Now, that you calculated descriptors for all the keypoinst, you have a way to compare those keypoints. For a simple example of image matching (when you know the images are of the same object, and would like to identify the parts in different images that depict the same part of the scene, or would like to identify the perspective change between two images), you would compare every keypoint descriptor of one image to every keypoint descriptor of the other image. As the descriptors are vectors of numbers, you can compare them with something as simple as Euclidian distance. There are some more complex distances that can be used as a similarity measure, of course. But, in the end, you would say that the keypoints whose descriptors have the smallest distance between them are matches, e.g. same "places" or "parts of objects" in different images.

For more complex use of keypoints/descriptors, you should take a look at this question -- especially the "low-level local approach" in my answer and "Bag-of-words" approach in @Maurits answer. Also, links provided in those answers are useful.


You show a pretty narrow notch, are you trying to notch out a tone? If so, maybe you could use an Inverse Chebyshev, and move the -10 dB pt in the notch to your tone frequency.


Properties of DTFT to Infer the Inverse DTFT of Altered Data

The property is the Linearity of the DTFT.

Linearity means that if your input is a linear combination of signals the output will be the same linear combination of each input by itself:

$$ \operatorname{DTFT} \left( \alpha x \left[ n \right] + \beta y \left[ n \right] \right) = \alpha \operatorname{DTFT} \left( x \left[ n \right] \right) + \beta \operatorname{DTFT} \left( y \left[ n \right] \right) $$

This also holds for the Inverse DTFT.

By the way, if you want people to keep answer your questions:

Mark them as solved once someone solves them for you.
Don't use images to post the questions. Write them so the question will be independent of out of site resources.


What is the best estimate of a signal given multiple noisy realisations

The term you're looking for is efficient estimator. That's the least-variance (== least RMS) estimator for a given amount of observation.
Whether or not that's the average depends on your noise distribution.
For example, consider the two noises:

continuous-valued uniformly distributed noise
Noise that only takes values $\in \{-\sigma, 0, \sigma\}$ with equal probability.

Both can have the same variance, both are zero mean. Both can be white. But while the average sounds like a good idea for the first one, for the second one you can get a 0 RMS estimator as soon as you've observed three different values by picking the middle one.


What is the difference between (Sample Timing Offset), (Carrier Frequency offset), and (Sampling Frequency Offset)?

I would have the easiest time explaining these briefly with the aid of an eye diagram as depicted in the graphic below.

In this example we see the constellation pattern of a raised cosine QPSK waveform on the complex plane on the left with I as the real axis and Q as the imaginary axis, and the resulting eye diagram pattern of the resulting real (I) and imaginary (Q) waveforms each on the right.  The eye diagram is what we would see if we triggered an oscilloscope to the symbol clock with infinite persistence, and here we see a 3 symbols duration superimposed over many symbols.

Notice the red dots on the left-- these are samples with no Sample Timing Offset, corresponding to the time locations on the right where the vertical variation between all the samples is minimized-- the ideal time to sample the waveform in order to make the final decisions on what data was transmitted. As we shift that sample clock to the right or left (Sample Time Offset), the vertical variation will increase, reducing our margin against noise to make correct decisions.

This constellation as shown is with no Carrier Frequency Offset (The waveform is centered at DC, the baseband signal). If there was a Carrier Frequency Offset, the entire constellation would rotate at the rate of that Offset-- In fact we measure the rate of rotation of the red dots to determine and correct for carrier offset error.

Finally Sampling Frequency Offset would be if the sampling clock itself had an error in frequency. This will change the effective Sample Time Offset and Carrier Frequency Offset proportional to the sampling rate error. If the sampling rate error was 1ppm, then the Sampling Frequency Offset and Carrier Frequency Offset as observed would also be in error by 1ppm. (For example, a 10MSps sampling clock with a 1 ppm error would have a 10 Hz error. If the observed carrier offset was 1 Hz with the 10MSps clock with no error, the carrier offset as observed would have an additional error of 1e-6 Hz.)

Note that in OFDM, a Sampling Time Offset can produce results similar to Sampling Frequency Offset in that the individual constellations in each sub-carrier will rotate from sub-carrier to sub-carrier (so a rotation in frequency rather than a rotation in time). This is because a time offset is a linear phase in the frequency domain. So pay very close attention to what domain you are in: If you see a rotation in the frequency domain, this is due to an offset in time. If you see a rotation in the time domain (on a single sub-carrier from sample to sample in time) then this is due to an offset in frequency.


This depends a bit on what you are after. As Laurent Duval has pointed out, you need can only take the logarithm of a unitless quantity so you need to do some form of normalization

If you want Sound Pressure Level (dbSPL) you would normalize to the standard reference pressure of 20 uPa (20e-6 Pa)
You could also divide by the free field acoustic impedance of air, which would give you the intensity in W/m^2.
The most common units for a PSD would be W/Hz. To calculate that you'd have to integrate over the entire volume flow and so you need to make an assumption for the shape of your sound wave. That's pretty easy for spherical and plane waves, but difficult for more complex shapes.


What is the significance of a spike like autocorrelation function

Consider a linear time-invariant system with unknown impulse response $h(t)$. If the input is a zero-mean wide-sense-stationary random process $\{X(t)\}$, then, with $\{Y(t)\}$ denoting the output process, we have that
$$R_{Y,X}(\tau) = E[Y(t)X(t+\tau)] = h\star R_x\big|_\tau.$$
So. if $\{X(t)\}$ is a white noise process wth autocorrelation function
$K\delta(t)$, then $R_{Y,X}(\tau) = Kh(t)$. Why is this at all important? Well, most processes are assumed to be ergodic, and so $R_{Y,X}(\tau)$ can be estimated by cross-correlating a sample path of the output with a sample path of the input with a delay of $\tau$. Thus, the spikiness of the autocorrelation function has a role in system identification of unknown systems.