For a given threshold $\lambda$ (that can be dependent on resolution level), and value of wavelet coefficient $d$, hard thresholding is defined as:

$  D^H(d|\lambda)=\begin{cases}
               0,& \text{for } |d| \leq \lambda\\
               d,& \text{for } |d| > \lambda
            \end{cases}$

whereas soft thresholding is governed by following equation:

$  D^S(d|\lambda)=\begin{cases}
               0,& \text{for } |d| \leq \lambda\\
               d-\lambda,& \text{for } d > \lambda \\
               d+\lambda,& \text{for } d < -\lambda \\
            \end{cases}$

Figure below depicts both cases:

The soft thresholding is also called wavelet shrinkage, as values for both positive and negative coefficients are being "shrinked" towards zero, in contrary to hard thresholding which either keeps or removes values of coefficients.

In case of image de-noising, you are not working strictly on "intensity values", but wavelet coefficients. You probably remember that you can decompose your image into wavelet levels, like in case of lovely Lena.

 Assuming that wavelet transform gives sparse coefficients, mostly close to zero, and noise level is lower than wavelet coefficients, you can simply threshold these. Although if you wish, you can perform hard/soft thresholding on each decomposition level with a different value of $\lambda$. When it is done, then you just have to reconstruct your image from all decomposition levels and voila, noise should be removed!

Below you have two examples of de-noised image via hard and soft thresholding respectively (same $\lambda$). Obviously soft thresholding gives more smooth image- if you can notice that with such a poor resolution ;) Courtesy of MATLAB.


Power Spectral Density of discrete White Gaussian Noise defined with variance and sampling interval

Marcus is correct, so this means something else is wrong. The left equality simply shows that, for continuous time, the unilateral PSD, $S_x(f)$, is two times the Fourier transform of the autocorrelation function, $R_x(\tau)$. See, e.g., A.B. Carlson, Introduction to Communication Systems, 2nd Ed., McGraw-Hill, NY ©1975, Chapter 2, equation 18a. Or see pretty much any other book on Fourier transforms.

In discrete time, with consecutive samples spaced by $\Delta t$, the unilateral PSD for white noise is $2\sigma^2\Delta t$. Small values of $\Delta t$ are usually desirable, since the sampling frequency is the reciprocal of $\Delta t$ and the Nyquist frequency is half the sampling frequency. But $\Delta t$ cannot go to zero, as Marcus commented.

What happens if you just delete the limit stuff, on the right hand side of equation 7.22, and see what happens in the textbook's subsequent matched filter treatment? For a matched filter, the maximum signal-to-noise ratio is simply twice the pulse energy divided by the unilateral (white noise) PSD. (Same reference as above, Chapter 4, equation 12a.) So it looks like $2\sigma^2\Delta t$ will be used in your textbook. Hope this helps a bit.


Matched filter for "amplitude SNR" vs power SNR

Power is proportional to the square of amplitude, so there's really no difference (since SNR is always a positive number) when maximizing: A filter that maximizes "amplitude SNR" maximizes "power SNR".

So, if your application requires you to maximize the amplitude ratio between signal and noise (be careful, amplitudes can be negative or complex, depending on how you define them!), then that's the same as requiring you to maximize the power ratio.

I.e.

$$h_\text{opt} = \arg\max_h E\left\{\left(h^H (vv^H)h\right)^{\frac12}\right\} = \arg\max_h E\left\{h^H (vv^H)h\right\} $$

since the quadratic form $h^H (vv^H)h$ (which is the power of $y$) can never be negative, and that makes your square root both unambiguously invertible and monotonous.


Uniqueness of Fourier Series Representation and the Fourier Transform of Periodic Signals

The Fourier series, if it exists, is unique, apart from the fact that we can always use a multiple of the fundamental period $T$ as the period, or, equivalently, a fraction of the fundamental frequency $\omega_0$. The only consequence of this is that we introduce Fourier coefficients with value zero in between the non-zero coefficients. The values of the non-zero coefficients remain unchanged:

$$x(t)=\sum_{k=-\infty}^{\infty}a_ke^{jk\omega_0}=\sum_{k=-\infty}^{\infty}\tilde{a}_ke^{jk\omega_0/N},\quad N\in\mathbb{Z}^+$$

$$\tilde{a}_k=\begin{cases}a_{k/N},\quad k=lN,\quad l\in\mathbb{Z}\\0,\quad\text{otherwise}\end{cases}$$


Using ROC curves for comparing the performance of SIFT and SURF

The ROC:
ROC curves are popularly used as performance metric for classification tasks. If the images in your dataset has class labels, then you can employ supervised learning to train a classifier (SVM for example). The dataset is split into training and testing and predicted class score from the classifier for images in the test set are compared to ground truth labels. The ROC curve is acquired by applying a threshold value to the classifier predicted score and obtaining a (TP,FP) value for each threshold to generate the curve.

The features:
SIFT and SURF are detectors and image local-patch descriptors (you can also use these without the detector in what is now more successful approach called 'dense sampling' where the keypoints are uniformly sampled in the image).

There are several approaches, but the most popular and recommended for first use is the 'Bag-of-Words' (BoW) model based approach. The reason is that: each image has one label (ignoring multi-class and multi-label scenario for now) and several hundred keypoints and each keypoint has an associated high-dimensional feature vector. The feature for an image at this stage is a matrix and we would like a vector (a point in feature space for each image with its label). The BoW approach does just this job for us.

BoW approach:
All the feature vectors for every image in the training set are collated and clustered (ignoring labels and location data) to compute a 'dictionary', using unsupervised learning. The dictionary is a set of cluster centroids, most commonly computed using k-means clustering. The dictionary follows data compression scheme to succinctly represent the data distribution.

Subsequently, each image is 'encoded' using its keypoints and the computed dictionary. Using vector quantization VQ, each keypoint is associated with a dictionary element to compute a occurrence histogram of keypoint associations. The consequence is that each image now has a feature vector which is the normalized occurrence histogram of keypoints assigned to dictionary elements. The dimensionality of this encoded feature vector is the size of the dictionary.

This normalized histogram for each image along with the label is used to training a classifier and then acquire classification results measured in terms of ROC curve. The sequence is:

[Images] -> (SIFT/SURF) -> [SIFT/SURF feature matrix] -> (collate and cluster) -> [dictionary] -> (BoW histogram by VQ encoding) -> [feature vector, label] -> Classifier training (training data) -> Classifier testing (validation/test data) -> ROC curve.

Suggestions:
Besides ROC, there are other performance measures available like accuracy, mean average precision, F1 score, etc. Your choice depends on the imageset, label ratio balance, objective of the classifier, etc.

There is plenty of code available online for baseline Bag-of-words model for images. One sample in Matlab using the Caltech101 dataset is at:
http://www.vlfeat.org/applications/apps.html#apps.caltech-101

The performance of SIFT and SURF is also provided in the original paper on SURF, so you can compare your results if you choose to work with the same dataset.


When I double bit rate I do not see a difference in I_spectrum and Q_spectrum, but why?

Do you know how ADS is plotting the spectrum? Plotting the spectrum without doing some kind of normalization will give you a higher magnitude. Once you determine the size of the FFT, normalizing by the length will give you the same magnitude no mater what sampling rate you choose.
For example let's take two rectangular signals, one sampled at 1 MHz and the other at 2 MHz. Below are their spectrums without normalization:

Since the bottom one is sampled twice as fast, it eventually produces an FFT size that is twice as long, hence the 6 dB increase in the peak.
Now compare this to the same exact signals, but now their magnitudes are normalized by their respective FFT sizes:

Now you can see that the peaks are the same magnitude. You can play with normalizing all day long to fit your need. It is the shape of the spectrum that is usually most important.
Here is some quick MATLAB code so you can maybe try it yourself and play around a bit.
%% Signal generation and FFT

% Sampling rates
fs1 = 1e6;
fs2 = 2e6;

% Rectangular pulse signals
t1 = 0:1/fs1:1e-5;
t2 = 0:1/fs2:1e-5;
pulseSignal1 = ones(1, numel(t1));
pulseSignal2 = ones(1, numel(t2));

% FFT setup
nfft1 = 100*numel(t1);
f1 = fs1.*(-nfft1/2:nfft1/2-1)/nfft1;

nfft2 = 100*numel(t2);
f2 = fs2.*(-nfft2/2:nfft2/2-1)/nfft2;

%% Without Normalization

figure;
subplot(2, 1, 1);
plot(f1./1e6, 20*log10(abs(fftshift(fft(pulseSignal1, nfft1)))));
xlabel("Frequency (MHz");
ylabel("Magnituide (dB)");
legend("F_s = 1 MHz");
ylim([-40 50]);

subplot(2, 1, 2);
plot(f2./1e6, 20*log10(abs(fftshift(fft(pulseSignal2, nfft2)))));
xlabel("Frequency (MHz");
ylabel("Magnituide (dB)");
legend("F_s = 2 MHz");
ylim([-40 50]);

%% With Normalization

figure;
subplot(2, 1, 1);
plot(f1./1e6, 20*log10(abs(fftshift(fft(pulseSignal1, nfft1)./nfft1))));
xlabel("Frequency (MHz");
ylabel("Magnituide (dB)");
legend("F_s = 1 MHz");
ylim([-80 -10]);

subplot(2, 1, 2);
plot(f2./1e6, 20*log10(abs(fftshift(fft(pulseSignal2, nfft2)./nfft2))));
xlabel("Frequency (MHz");
ylabel("Magnituide (dB)");
legend("F_s = 2 MHz");
ylim([-80 -10]);


Power spectral density of sigma-delta modulator

Everything you have written sounds correct so I assume you forgot to consider the BW for each FFT bin and are comparing to an expected spectral density that is per Hz.  Without windowing or zero padding the BW in each bin will be 1/T where T is the time length of your signal (so in other words 1 bin wide).  Therefore you want to reduce the total power measured in each bin by 10Log(1/T) to get the equivalent power per Hz spectral density.

If you are using windows you can refer to fred harris' classis paper on windowing for the Equivalent Noise BW of all common bins as well as windowing implementation loss for an accurate computation.

https://www.utdallas.edu/~cpb021000/EE%204361/Great%20DSP%20Papers/Harris%20on%20Windows.pdf


Multi Image (Same Noise Realization) Denoising

For me, if all you want to deal with is the possible 'debris in the optics,' why not get a reference background image? This will help you to easily identify which are 'noise' caused by your optics. Assume you can get a perfect background, namely everything remains the same except for your cells ( maybe you can get this by taking another image after removing all cells ), then simply subtract the background image from the original one will work.


Just to follow up on this, I believe I figured out one way to do what I was looking for shown below...if anyone has any other ideas feel free to share too


Which method should be selected to enhance the compressed image?

I would advise you to try between two different approaches, then test and see which one performs best.

The first one is intensity clustering. In this approach, you data points are one-dimensional and correspond to the intensity of your raw pixels. Then, run a clustering algorithm such as $k$-means with $k = 256$, and perform the final quantization by replacing the intensity of each pixel by the mean value (or the index) of the corresponding cluster mean. With this approach you can limit the information loss by forming an image of indices and keeping a look-up table of the real (14 bits) cluster mean values.

The second would be to perform some local histogram equalization, i.e., cutting your image into blocks and performing equalization (or any other enhancement procedure) on each block.
You will obtain an image where the features are mstly preserved and that you can remap linearly into $[0,255]$.


If compute time doesn't matter, you can overcome the memory problem easily with a dual for-loop implementation. But there's an FFT alternative with caveats: break up signal into M seqments, take M FFT's, add them up. The result is exact, limited in the lowest possible frequency it can represent. You could, however, do both; FFT for high frequencies, and for-loop for leftovers.

FFT in segments
Recall that FFT computes DFT, which in turn multiplies input by sines of cosines of varying frequencies. Let's take x of length 512, f=4 for left half, f=64 for right:

real(FFT(x)[1]), for example, was computed by multiplying:

And imag(FFT(x)[1]) would be computed by replacing the orange cosine of f=1 with orange sine of f=1. Now, what if N=512 is too long, and we decide to break it up into M=2 frames? Then, from those FFT's (FFT256) perspective, real(FFT(x)[1]) looks like

and from the FFT512's perspective, this looks like

which is real(FFT(x)[2]). In other words, FFT(x[:256])[1] + FFT(x[256:])[1] == FFT(x)[2].
Now imagine the highest frequency basis (sine or cosine) in either FFT256 or FFT512; what's the limit? It's [-1, 1, -1, 1], sample-to-sample (k=1/N). No matter how long the FFT is, this basis is always the greatest (normalized) frequency possible from x's point of view. As from above example, shortening segments corresponded to being unable to compute FFT(x)[1] - because that'd require a k=0.5 basis in FFT256:

It's important to understand "normalized frequency" to correctly combine the FFT's later. But basically, numerically the shorter segment computations are exactly the same as longer segment computations with higher-frequency bases.
For a billion point signal, the highest FFT basis frequency is k=500,000,000. What this actually means, is that the dynamic range of frequencies you can capture is 500,000,000:1 (look up normalized frequency). So, if in data you expect the ratio between the highest and lowest frequency to be 10 times less, then you can capture the entire spectrum by splitting up the billion points into 10 segments and adding up the results. If it's any more, you can capture a part of it with the above method, and another part with for-loops.

Additional notes:

If the input is real-valued, be sure to use rfft (real FFT), which takes half the memory.
Beware of any solutions proposing a window of any kind other than rectangular; they will distort the spectrum (but it may be practically insignificant).
Beware zero-padding
For very long DFT, be sure to use float64, and even that may not be enough depending on level of accuracy you seek

For-loops rDFT implementation: (Python but easy to port to MATLAB)
def dft(x):  # unnormalized
    N = len(x)
    reX = np.zeros(N // 2 + 1)
    imX = np.zeros(N // 2 + 1)

    for k in range(N // 2 + 1):
        for n in range(N):
            reX[k] += x[n] * np.cos(2 * np.pi * k * n / N)
            imX[k] += x[n] * np.sin(2 * np.pi * k * n / N)
    return reX - imX * 1j

Additionally
If working in Python, I recommend Numba for CPU-parallelism, and Cupy to run code on GPU. Both help with overcoming Python's slug for-loops and Numpy's single-threadedness.


Acceptable Carrier Frequency Offset for Specific BPSK system

Bottom line: Either implement a (simple) carrier recovery algorithm or use DBPSK instead of BPSK if you are ok with the 3 dB SNR penalty in performance.

If you want to do synchronous detection of BPSK, you will need to completely remove the carrier offset. If you do not correct for carrier offset, then your BPSK symbols will continue to rotate, eventually inverting; unless you know the rate of inversion you will not be able to distinguish if the data transmitted was a 1 or a 0. Knowing the rate of inversion is essentially carrier recovery/tracking.

See the figure below specific to your case showing the transmitted data in blue at your 1 KHz data rate, and the resulting real portion of your received signal with a 60 Hz carrier offset.

Even if you reduced this carrier offset to significantly small values, you will still reach a point in time where the data that was completely real moves to the imaginary axis and then becomes real again but will be inverted. In the simplest case to answer your question, if your frequency offset was sufficiently smaller than the length of your data packet (768 bits), then it is feasible that you could determine the inversion from a preamble (but you could still have an inversion in the middle of your data depending on your relative phase, but let's assume you do a simple phase correction in your preamble as well; now if you do both phase and frequency correction in your preamble- then you are doing carrier recovery). So in this case the maximum offset would be:

$$f_{offset} < \frac{1000}{2\times768} = 0.65 \text{ Hz} $$

Which is the frequency above which the frequency offset would have rotated your data more than 180 degrees over the duration of the packet.

You can certainly implement additional coding in your message payload and therefore allow for even greater offsets-- however DBPSK would be a much easier solution on this front, but coding is feasible if you knew you could guarantee clock accuracy and wanted better SNR performance at the expense of some overhead). However be aware that with a separate transmitter and receiver clock that even if you were to manually adjust them so that you knew the clock offsets were below these limits, those clocks will drift over time and temperature where you will not be able to guarantee with your rates that the clocks will remain synchronized open loop. For example, the reasonable drift rate between two separate 10 MHz crystal oscillators is 100 Hz/day. With wide temperature changes this can be significantly worst.  That said, you could periodically measure your frequency error and from knowledge of the worst case frequency drift for your design, this measure (and correction) rate could be minimized for maximum efficiency. This is still a carrier recovery implementation with some suggestions detailed below on that front, but is to suggest that the carrier recovery operation could feasibly be once per packet, or even significantly less when you are resource constrained.

DBPSK

DBPSK (Differential BPSK) is an non-coherent BPSK modulation, meaning we do not require the transmitter and receiver to be synchronized. The block diagram showing the DBPSK demodulation process is shown in the figure below. T is the delay of one symbol, so we compare the phase of one symbol to the previous symbol to determine if a data transition occurred (a multiplier followed by a low pass filter is a phase detector with the response $V_{out} = Acos(\theta)$, so if two adjacent symbols are in phase (at any phase angle!) the output will be +A, and if they are 180° out of phase, the output will be -A. If you were not concerned with establishing the data boundaries (meaning something else downstream is doing this for you), then you could use a simple low pass filter (a moving average over a symbol duration T is a great choice due to its simplicity and the fact that it will have frequency nulls at 1/T), otherwise you will need to also establish the locations of the symbol boundaries (timing recovery) and then an integrate and dump would be the optimum solution for data decision.

Differential encoding https://en.wikipedia.org/wiki/Differential_coding is often used in conjunction with DBPSK since the demodulation is only telling us if a symbol changed or not, but will have arbitrary inversion depending on which symbol was processed first. The demodulator is in fact providing the decoding function of differential encoding.  So in such a system, as long as the relative clocks of the transmitter and receiver are not drifting to cause more than a 180° over the time duration of 2 symbols, there will be no frequency offset related issue with demodulation. This means in your example, that over 2 ms the frequency of the transmitter or receiver (and they can be completely different frequencies) cannot change to create an additional phase offset of 180 degrees or more, which leads to a change in frequency of $\frac{1}{0.004} = 250 \text{ Hz}$. So as long as both your clocks are stable such that they do not change more than 250 Hz in 4 ms you will have no issue with this approach (suggesting that you could even use free-running undisciplined local oscillators for an ultimate low cost approach).  The primary drawback of this approach is you are using your own signal as the local oscillator, so any noise in your signal gets added to itself in the process (doubling the amount of noise for any noise that is uncorrelated from symbol to symbol). That is why DBPSK has a 3 dB SNR penalty.

Carrier Recovery

See this post for common implementations of carrier recovery for BPSK: Phase synchronization in BPSK and I draw your attention to this phase detector implementation specifically that is shown in the figure below:

If in your implementation you have the complex I and Q samples for your received data, then this could lead to a very simple carrier detection and removal technique since it only requires two real multiplications. Further, you could minimize how often you need to make this measurement based on your time drift and SNR constraints (if not resource limited, computing on every sample and averaging via the recovery loop would be optimum).

Otherwise other common techniques such as squaring your signal (which will give you a unmodulated tone that is exactly twice your carrier), and then dividing that by two to create a local synchronized carrier can be done which is also described in the referenced post.


When can the $\mathcal Z$-transform be inverted? When not?

If I understand your question correctly then you're asking under which conditions the inverse $\mathcal{Z}$-transform of a given function $F(z)$ exists. Since the inverse $\mathcal{Z}$-transform is just a Laurent expansion of $F(z)$ around $z=0$, the inverse transform exists if $F(z)$ can be developed in a Laurent series about $z=0$. This is the case if $F(z)$ is analytic in an annulus centered at $z=0$.


Modeling multipath distortion and shifting by fractions of a sample

One approach is to upsample the signal, then delay it, then downsample (or, if feasible, you can do the rest of the processing at the higher sampling rate).

This general approach to modeling multipath is correct. However, you need to make sure that your model is realistic and/or useful for your purpose. A few things to consider:

If you add delayed/attenuated replicas of the signal to itself, then you're assuming you have line-of-sight between the transmit and receive antennas. The effect is different when you don't have line-of-sight.
The number of paths is important too. A channel with many paths is easier to understand theoretically, because you can make certain assumptions about the overall channel gain (e.g. that it's a Gaussian random variable).
The maximum delay is important too. A very large delay introduces very destructive effects on the received signal. Many, very small delays can be modeled as a single gain factor.
In a static environment, the channel is also static. If the antennas and/or obstacles in the environment are in movement, then the channel is dynamic, in the sense that the number of paths, and the gain of each, vary with time.


Signal processing becomes more valuable and useful when applied into the problems of different fields such as mining, seismology, oceaonography, climatology etc (apart from the obvious applications in multimedia, communications and control).

There are devices used to test and measure strength, stress and fracture behaviour of materials that's to be encountered in nature or in man made structures so that a good understanding of the materials stress strain behaviour enables stronger, durable and safer constructions.

One method of obtaining such important information of materials comes from specific lab equipments that test a sample material by providing a controlled force and measure resulting displacements. And the material characteristics is analysed according to the recorded force and displacement measurements.

Therefore the method requires these measurements to reflect the true conditions as accurate as possible in order to make reliable and precise conclusions from them.

There are many reasons why the measurements can be imperfect however. The distortions must have been analysed considering the sample-device mechanical dynamic interactions and any electronic noise or device nonlinearities that can also interfere with the measurements.

Under these circumstances, consider the following imperfect measurement results of a given sample test run:

As can be seen, this is quite an unexpected result from seemingly such a stable (slowly moving) device. One would for example expect a more monotonic curve (though not necessarily a strictly monotonic curve but at least not so much oscillatory either...)

Looking in detail to the first 8000 samples gives the following plot:

A Fourier analysis at this point could be of great help, revealed by the following plot:

Now we are enlightened by the presence of those harmonic spikes that are present at about 8 Hz, 16 Hz and 24 Hz respectively. In addition, there is also a DC looking spike at the origin, but that's actually the signal of interest that we consider as the true measurment which is very slowly (hence DC looking) moving down from -0.5 mV to -1.5 mV on the previous figure.

To verify the presence of 8, 16 and 24 Hz oscillations we further zoom into the signal to observe them directly in time domain rather than indirectly provided by the Fourier analysis computation...

In this very clear cropped portion of length about 100 samples, from the sample 1460 to sample 1560, there are approximately 8 periods of a fundamental oscillation. (And it's also understood from its waveshape that this is not a pure sine wave but involves harmonics of it which create the unsymmetric up down appearence).

Given the sampling period of $Ts=0.01$ seconds approximateley, this makes the fundamental period of oscillations to be $T_0 = 0.01*100/8 = 0.125$ seconds and therefore the fundamental frequency of the oscillations is found about $F_0 = 1/T_0 = 1/0.125 = 8$ Hz, Just as predicted from the Fourier analysis plot.

So therefore we are convinced both directly from the time domain measurement and also indirectly from the frequency domain computation that there are three quasi-oscillatory components at about 8Hz, 16Hz and 24Hz in addition to the signal itself around DC frequencies.

So in this case most probably due to mechanical interaction of the force piston and the sample material, there are uncontrolled vibrations occuring, which are casuing the oscillatory measurements. Or it could as well be electronic noise picked up by some reasons undecidable from the measurements alone.

Now, the are a number of solutions which would enable us to get rid of the unwanted oscillations (provied that these oscillations are indeed unwanted !) One of the simplest of them is to use a lowpass filter with an appropriate cutoff frequency. In the following I'm using a simple FIR lowpass filter of cutoff frequency below 8 Hz so that those oscillaitons are removed, but the useful signal is retained. The resulting plot is the following:

Clearly (or not!) this is an improvement, but a number of extra trials must have been performed in order to decide on the best solutions. The number of taps used was quite high approaching 256. This would induce a group delay of 128 samples, hence care must be taken to interpret the exact timing of the resulting samples. But this is not at all a problem, as long as you can perform non-real time (offline , non-causal) filering and that you know precisely how much delay there will be in the output which is the case here.

A much shother filter could have been designed which would yield much less sample shifts. Also instead of a lowpass filter a harmonic notch filter could also be designed, the decision of which is guided by the useful signal charactheristics.


How to Get Rid of Ripples from a Gradient Image of a Smoothed Image?

I think it happens due to 2 things:

Quantization
You are working using UINT8 Image, try convert it into floating Point Image.
You may do this by mO = im2double(mI) where mI is the UINT8 and mO is a floating point image in the range [0, 1]. You may also do it using mO = double(mI) / 255.
Gaussian Kernel Radius vs STD Ratio
Your Gaussian Filter Radius is 12 while $ \sigma = 4 $. This is a ratio of 3. Try higher ratios, something like kernelRadius = ceil(5 * kernelStd) this will make the Gaussian Filter result much smoother.


How to Solve Image Denoising with Total Variation Prior Using ADMM?

I have done a bit of this myself and you'd need to adapt.
There is a Douglas Rachford self implemented and a primal dual approach here implemented in Recovery of Fusion Frame Structured Signal via Compressed Sensing.
Note that Clarice Poon (Bath University) had some nice tutorials on it.
Another source is the Numerical Tours from Gabriel Peyre. See Denoising by Sobolev and Total Variation Regularization.


How should I remove obvious anomalies from spectral data?

This seems to be a textbook case for a narrow median filter. Your $S$ is very broadband compared to your $N$, so a median filter will not affect your signal, as it can be designed quite narrow. And by using median instead of averaging (i.e. lowpass) you minimize the error introduced by the anomalies. I would give it a try.


Not really. Welch's method is specifically designed to estimate spectral density and discards phase information.
IFT and FT are almost identical and have the exact same complexity. If you want to do an IFT, just do an IFT unless you have a really good reason not to.


What are the definitions of symbols and samples related to PCM?

Quantization and encoding are largely independent. "Symbols" is another word for "pulses", and the line encoding can also play a role in how the information is transmitted.

Say you quantize one sample to $256 = 2^8$ levels, or 8 bits/sample. In order to transmit those 8 bits, you can, among other options:

use binary encoding, which requires transmitting eight symbols, each with amplitude in  the set $\lbrace +, - \rbrace$, where $+$ and $-$ correspond to some predefined voltage levels.
Use ternary encoding with AMI (alternate mark inversion), which also requires eight symbols, but now with amplitudes in the set $\lbrace +, -, 0 \rbrace$.
Use quaternary (or 2B1Q) encoding, which requires four symbols.

Other examples may include Manchester encoding, CMI, block line codes.... but in all cases, one symbol correponds to one transmitted pulse.


Processing an ECG signal with a median filter

You are certainly doing something very wrong. You should upload your data in order to get better responses. You can upload your data into any upload sites you wish and provide its link here, practically.

Median filter is a highly highly nonlinear filter (it re-orders the sample positions!). The output of the median filter at a position $n$ is the median of the values that reside in the window scope; i.e., it's the value that resides in the middle when the samples are sorted in order. Hence median filtering requires sorting for each computation. This makes it quite slow as well (a deailed answer actually depends on the architecture...)

Median filter is mainly used for speckle or salt and pepper noise removal, in essence these are local noise samples whose frequency domain filtering is not possible without degrading the whole signal. Such local (in time) peaks will have wide band frequency spectrums which inhibit frequency domain attacks to remove them, therefore, leaving only the time domain (or time-frequency domain) approaches possible.

Median filter has a tendency to preserve edges, therefore quite preferred in certain image enhancement operations. However it also has the side effect of washed out results (texture details are lost and only strong edges remain) which indicates that their use should be performed with care.

In principle, the longer the window size, the stronger will be the washed out effect. So it's customary to use as short as possible window sizes (unless otherwise dictated by the particular application)

Coming to your example plot. Using a window size larger than 50 samples will wash out local details that might be important to you, so you should use a window size less than 50, I guess (based on the plot you provided).


Why is the last value of an RFFT always real?

The real component of an FFT or DFT represents cosine or even sinusoids (symmetric around the middle), and the imaginary component represents sine, or odd sinusoids (anti-symmetric around the center).  An anti-symmetric waveform has to be zero at the center of an even length window.  If a sinusoid of frequency Fs/2 is zero at any sample, it is zero at all samples taken at a sample rate of Fs.  The imaginary component of the basis vector at N/2 for a DFT transform is also all zeros for the same reason. Thus the imaginary bin at Fs/2 is the sum of zero times zero, and has to be zero.  Leaving the result in bin N/2 strictly real for even length DFTs.  Perhaps this is the bin your RFFT implementation reports for A(-1).

The is also the reason why Nyquist sampling has to be at a rate ABOVE the_highest_frequency / 2, not equal to.  Sampling will leave out part of any anti-symmetric waveform or waveform component (in any even length set of samples) if that component is of exactly frequency F when sampled at F/2.

There are different rules for odd length FFTs.


Difference between librosa.stft and scipy.signal.stft

Fft scaling (normalization) isn't strictly a part of an fft algorithm although some implementations (like the one you mention) include it.
Since scaling depends on the actual fft algorithm, it makes sense, in my opinion,  to include it in the implementation because it may not be obvious for programmers to decide what scaling factor(s) to use (the usual scaling factors are $1/N$, $2/N$ (an efficient fft algorithm that is using only half the spectrum so we need to multiply by 2 to compensate for the missing energy), $1/sum(wnd)$, $2/sum(wnd)$
As for your second question, I haven't gone over the code (you haven't provided any links to the source), but some (usually minor) calculation differences might exist
given that an fft algorithm can be implemented in a number of ways ( e.g real fft vs complex fft, decimation in frequency vs decimation in time,
radix 2 vs mixed radix etc and combinations thereof)


Does it mean that if the source transmits information to relay, it needs a carrier to help it send the message, so if I have K source, then I have K subcarrier?

There's $K$ pairs of nodes, if you read the text closely. That means you can have, for example, one source with four information sinks, or four sources and three sinks, or…


Signals are (almost) always time variant. Otherwise they would be very boring :-). I believe your question is about whether "systems" are time invariant or not. A system basically is a "prescription" of how to create an output signal from an input signal. The equation you show is a system equation of how to create output y[n] from input x[n].

The idea of time invariance is the following: The behavior of the system doesn't depend on absolute time, i.e. the relationship between input and output is the same, regardless of whether it's Tuesday or Friday. You can find a more mathematical description here http://en.wikipedia.org/wiki/Time-invariant_system

The easiest invariance test is the following:

pick any input signal, typically $\delta [n]$ will do fine.
Calculate the output $y{_{0}}[n]$ for that intput
now delay the input signal by some amount.
Calculate the output $y{_{1}}[n]$ for the delayed input
If $y{_{1}}[n]$ is NOT simply a delayed version of  $y{_{0}}[n]$ then the system time variant

Another (equivalent) test is the following

For time invariant system, the system equations only dependency on time must be as linear index into the input signal

Either one will work.


Good Language for Implementing Computer Vision System / Algorithm / Application?

Maybe you can be more specific about the scope and scale of your work (academic project? Desktop or Mobile commercial product? Web-based commercial project?).

Some recommendations and comments:

Matlab is common in the academic world, and quite good for sketching/validating ideas. You will have access to a large body of code from other researchers (in CV and machine learning); prototyping and debugging will be very fast and easy, but whatever you will have developed in this environment will be hard to put in production. Depending on what your code is doing, you might have memory/performance problems (there are situations where you can't describe what you want to do in terms of Matlab's primitives and have to start looping on pixels and Matlab's being an interpreted language is not helping in this context). Interaction with databases, web servers etc is not easy, sometimes impossible (you won't get a Matlab program to become a Thrift server called by a web front-end). Costs $$$.
C++ is what is used for many production-grade CV systems (think of something at the scale of Google's image search or Streetview, or many commercial robotics applications). Good libraries like OpenCV, excellent performance, easy to put into a production environment. If you need to do machine learning, there are many libraries out there (LibSVM / SVMlight, Torch). If you have to resort to "loop on all pixels" code it will perform well. Easy to use for coding the systems/storage layers needed in a large scale retrieval system (eg: a very large on-disk hash map for storing an inverted index mapping feature hashes to images). Things like Thrift / Message Pack can turn your retrieval program into a RPC server which can be called by a web front-end. However: not very agile for prototyping, quite terrible for trying out new ideas, slower development time; and put in the hands of inexperienced coders might have hard to track performances and/or instability problems.
Python is somehow a middle ground between both. You can use it for Matlab style numerical computing (with numpy and scipy) + have bindings to libraries like OpenCV. You can do systems / data structure stuff with it and get acceptable performances. There are quite a few machine learning packages out there though less than in Matlab or C++. Unless you have to resort to "loop on all pixels" code, you will be able to code pretty much everything you could have done with C++ with a 1:1.5 to 1:3 ratio of performance and 2:1 to 10:1 ratio of source code size (debatable). But depending on the success of your project there will be a point where performance will be an issue and when rewriting to C++ won't be an option.


"Low-pass filter" in non-EE, software API contexts

The filter in your example is a first-order infinite impulse response (IIR) filter. Its transfer function is:

$$
H(z) = \frac{1 - \alpha}{1 - \alpha z^{-1}}
$$

which corresponds to a difference equation of:

$$
y[n] = \alpha y[n-1] + (1-\alpha) x[n]
$$

where $x[n]$ is the filter input and $y[n]$ is the filter output.

This type of filter is often used as a low-complexity lowpass filter and is often called a leaky integrator. It is favored because of its simple implementation, low computational complexity, and its tunability: its cutoff frequency depends upon the value of $\alpha$. $\alpha$ can take on values on the interval $[0,1)$. $\alpha = 0$ yields no filtering at all (the output is equal to the input); as $\alpha$ increases, the cutoff frequency of the filter decreases. You can think of $\alpha = 1$ as a boundary case where the cutoff frequency is infinitely low (the filter output is zero for all time).

You can think of this intuitively by noticing that the filter input is weighted by $\alpha$, so as the parameter increases, the quantity $1-\alpha$ decreases, so each input sample has a smaller proportional effect on the value of any particular output sample. This has the effect of smearing out the filter's impulse response over a longer period of time. Summing over a longer period of time is similar to computing a long moving average. As the length of a moving average increases, the cutoff frequency of the average decreases.

For your example, where $\alpha = 0.8$, the frequency response of the filter is as follows:

From the example, I would guess that this filter is being used to smooth high-frequency noise out of a time series of measurements from a sensor, trying to tease out a comparatively low-frequency signal of interest. This would be a very typical application for this sort of filter.

On your other sub-question, you are correct that filtering is often implemented via convolution of the input signal with the filter's impulse response. In most cases, this is only done with finite impulse response (FIR) filters. IIR filters such as this one are typically implemented using the filter's difference equation; since an IIR system's impulse response is infinitely long, you must truncate it to some finite length to make convolution with it tractable, at which point the filter is no longer IIR. The difference equation format is almost always cheaper to implement computationally, although the feedback inherent in that structure can lead to numerical issues that must be addressed (such as internal overflow and roundoff error accumulation).


Calculating the phase shift between two signals based on samples

It's an approximation as also admitted by your statement "satisfactory performance". Here is why.

Let $a[n]=\cos(w_0 n)$ and $b[n]=\cos(w_0 n + \theta)$ be two sequences of size $N$, represented by $1 \times N$ row matrices $\bar{a}$ and $\bar{b}$ in a program (such as Octave or Matlab)

Treating those two sequences $a[n]$ and $b[n]$ which are represented by two row matrices $\bar{a}$ and $\bar{b}$ as components of vectors $\vec{a}$ and $\vec{b}$, we can take advantage of the dot product stated for vectors and (approximately) compute the phase angle $\theta$ in between those two sinusoidal sequences $a[n]$ and $b[n]$

Geometric evaluation of the dot product between two vectors of same size is:
$$ \vec{a} \cdot \vec{b} = |\vec{a}||\vec{b}| \cos(\phi)$$

And the angle between those two vector therefore is:
$$\cos(\phi) = \frac{\vec{a}\cdot \vec{b}}{|\vec{a}||\vec{b}|}$$
Where $\phi$ is the angle between those vectors,  which will be equal to $\theta$ as shown:

The dot product between the vectors $\vec{a}$ and  $\vec{b}$ can also be algebrically computed from their MAC sum of components which are contained in the elements of the matrices  $\bar{a}$ and $\bar{b}$:
$$ \vec{a} \cdot \vec{b} = \sum_{i=1}^{N} \bar{a}(i)\bar{b}(i) $$

The product $\bar{a}(i)\bar{b}(i)$ can be shown to be equal to the following (using trigonometry):
$$\bar{a}(i)\bar{b}(i) = \cos(w_0 i) \cos(w_0 i + \theta)= 0.5\cos(2 w_0 i + \theta) + 0.5 \cos(\theta)$$

And the dot-product summation becomes:
$$ \vec{a} \cdot \vec{a} = \sum_{i=1}^{N} a(i)b(i)  =  0.5 N \cos(\theta) + \sum_{i=1}^{N} 0.5 \cos(2 w_0 i + \theta)$$

The absolute values $|\vec{a}|$ and $|\vec{b}|$ of the vectors $\vec{a}$ and $\vec{b}$ can also be computed from their norms computed over the matrices $\bar{a}$ and $\bar{b}$ as follows:

$$ |\vec{a}| = ||\bar{a}||=  ( \sum_{i=1}^{N} a(i)^2 )^{0.5}$$
$$ |\vec{b}| = ||\bar{b}||=  ( \sum_{i=1}^{N} b(i)^2 )^{0.5}$$

Again expanding the squared terms result in:
$$ |\vec{a}| = ||\bar{a}||=  ( 0.5 \sum_{i=1}^{N} 1 + \cos(2 w_0 i)  )^{0.5}$$
$$ |\vec{b}| = ||\bar{b}||=  ( 0.5 \sum_{i=1}^{N} 1 + \cos(2 w_0 i + \theta) )^{0.5}$$

Now, it can easily be shown that the following summations are small compared to large $N$:
$$ N+\sum_{i=1}^{N} \cos(2 w_0 i) \approx N , ~~~ N+\sum_{i=1}^{N} \cos(2 w_0 i + \theta) \approx N  $$
for large $N$ and suitable $w_0$:

Simplify the above equations to get the angle:
$$\cos(\phi) = \frac{\vec{a}\cdot \vec{b}}{|\vec{a}||\vec{b}|}$$
$$\cos(\phi) \approx \frac{ 0.5 N \cos(\theta)}{|\vec{a}||\vec{b}|}$$
$$\cos(\phi) \approx \frac{ 0.5 N \cos(\theta)}{0.5 N}$$
$$\cos(\phi) \approx \cos(\theta)$$
$$\phi \approx \theta$$

The approximation depends on the length $N$.


Choosing the best measure for machine learning on a Spectrogram

First thing: a spectrogram does not give the "frequencies" at time $t$. It gives amplitudes, measured for a set of linearly spaced frequency bands.

Anyway, there is no "right" answer to your question, because it depends on which musicological aspect you want to focus on, and you need to enlighten us on that. If you want to focus on melody, extracting the fundamental frequency is the most relevant feature, but I suspect this is not the case here. If you want to discriminate different timbres, spectral centroid might be a better way of "summarizing" a slice of spectrogram into a single value - better, but still quite bad.

Indeed I can see several flaws in your methodology:

Why extract information at time t only (assumedly at a small window centered at t) rather than on the entire segment of sound between onsets?
Why do you have this constraint that for any onset you have to extract a single real value? Some properties of sound, like timbre, are multidimensional, so you won't be able to entirely capture them in one single value. If we take the example of discriminating tabla sounds (which might or might not relevant to your problem, since you mentioned tala...), we need several features to discriminate them (bayan/dayan discrimination relies on pitch information, while resonant vs non-resonant stroke discrimination relies on timbral features).
You are completely discarding the rhythmic (time) information. Is that really wanted for your problem?

Tell us more about your problem, there's a good chance it has already been solved, and that you are asking the wrong questions to solve it...


Ambiguous results of phase correlation in pattern recognition

I suspect your problem occurs due to some scaling issues. Basically you need to normalize your research image to the pattern template by subtracting the mean value of the template. And it is better calculate the ratio of correlation to the standard deviation of both images.

I don't know which programming language you are using. I wrote a Matlab code for you to better get what I meant above (note that I did not use zero padding since I work on the spatial domain, I don't think it is the issue that caused your spurious correlation too):

file1='https://lh3.googleusercontent.com/-GQSrln7qnO8/UtCT8elB_FI/AAAAAAAAEaw/bh7aKyyMRIo/s512/subimage.png'; % Research image
I1=double(rgb2gray(imread(file1)));
file2='https://lh5.googleusercontent.com/-YaDK_1QvqIQ/UtCT6-rNkNI/AAAAAAAAEao/41MOQJ_bXxU/s512/pattern_grey.png';   % Template pattern image
I2=double(rgb2gray(imread(file2)));

I2=I2(:,1:32); %% No padding

It=double(I2);  % template
Ii=double(I1);  % image

Ii_mean = conv2(Ii,ones(size(It))./numel(It),'same');
It_mean = mean(It(:));
corr_1 = conv2(Ii,rot90(It-It_mean,2),'same')./numel(It);
corr_2 = Ii_mean.*sum(It(:)-It_mean);
conv_std = sqrt(conv2(Ii.^2,ones(size(It))./numel(It),'same')-Ii_mean.^2);
It_std = std(It(:));
S = (corr_1-corr_2)./(conv_std.*It_std);

imagesc(abs(S))

What you get:

You can observe the correlation peak at 17th row, 450th column, which is the position of your pattern in original image.

You can also use Matlab built-in function normxcorr2(It,Ii);, just note that this function results in the correlation map with the sum of the lengths of two images, while my code provides you the correlation map the same size of your research image.


Why is this recursive DFT algorithm not equivalent to this iterative method?

All right, you made two small mistakes and here is the corrected version: (I changed your function's name). Btw, you can get rid of the for loop too... ;-)

function [ s ] = recDFT2( x )
%
% a recursive DFT
%

N = length(x);

if N == 1
    s = x(1);
    return;
else

    X1 = recDFT2(x(1:2:N));
    X2 = recDFT2(x(2:2:N));
    s = zeros(N,1);

    for k = 1:(N/2)
        W(k)     = exp(-1j*2*pi*(k - 1)/N);
        s(k)     = X1(k) + W(k) .* X2(k);
        s(k+N/2) = X1(k) - W(k) .* X2(k);
    end
end


Why does the standard waveform from audio out get distorted when I reduce the amplitude?

Primarily because the documentation is poorly written. If you want to play at half amplitude you need to use

soundsc(y,48000,[-2 2]);

The range argument is the input amplitude range that gets mapped to the maximum available range of the sound card. So if the max of your input exceeds the max in the range argument, you will get clipping.


Well known waveforms or signal equations for filter testing

My question is whether there is a set of well known waveforms with algorithms that can be programmatically generated?

Yes and no. There are many wave forms that can be created programatically. However these may or may not be useful for you test.

Your test should focus on determining whether the algorithm meets the stated requirements. Hence the best way to build a test is to run down the list of specifications and requirements and build test that clearly establish whether the implementation meets the requirement or not . Once you know what exactly a specific test needs to assess and what the pass/fair criteria are, you can start thinking about how to build the test. That includes what wave forms but also what test equipment, what test flow, what data acquisition, what metrics etc.

Wave forms are a means to an end. Once you have a clear definition of what the end is, the wave form question will most likely become obvious.


No, Equation (2) is not the sum of 2 FFTs (of lengths $N/2$) as you claim, and it is not computationally more efficient. The canonical way of evaluating Equation (1) is via Horner's method:

$$X[k] = (((\cdots (x[N-1]\cdot\omega + x[N-2])\cdot\omega + x[N-3])\cdot\omega + \ldots)\cdot\omega + x[0]\tag{a}$$

where $\omega = W_N^k$ which requires $N-1$ multiplications and $N$ additions: lather, rinse, and repeat for for other values of $k$. Equation (2) just splits this Horner's rule computation into two parts with no change in complexity. That is, we compute $\sum_{n\ even}x[n]W_{N}^{nk}$  as $$(((\cdots (x[N-1]\cdot\Omega + x[N-3])\cdot\Omega + x[N-5])\cdot\Omega + \ldots)\cdot\Omega + x[0]\tag{b}$$ where $\Omega = \omega^2$ and
$\sum_{n\ odd}x[n]W_{N}^{nk}$ as
$$(((\cdots (x[N-2]\cdot\Omega + x[N-4])\cdot\Omega + x[N-6])\cdot\Omega + \ldots)\cdot\Omega + x[1])\cdot\omega.\tag{c}$$ Each computation requires $N/2$ multiplications and $N/2$ additions with no savings in computation over the standard brute-force Horner's method applied to a vector of length $N$ in $(a)$.

Equation (3), however, is the core of the FFT butterfly; the two sums in Equation (3) are DFTs of two shorter length sequences.  Define \begin{align}x_{\text{even}}[k] &= x[2k], &k = 0, 1, \ldots, N/2-1,\\
x_{\text{odd}}[k] &= x[2k+1], &k = 0, 1, \ldots, N/2-1,
\end{align} and let the DFTs of these shorter sequences be denoted by $X_{\text{even}}[k]$ and $X_{\text{odd}}[k]$ respectively. Then the computation of Equation (1) has been re-written in Equation (3) as
$$X[k] = X_{\text{even}}[k] + W_N^kX_{\text{odd}}[k] \tag{**}$$
which requires one multiplication and one addition per bin. Note that $k$ varies from $0$ to $N-1$ in $(**)$, and the shorter DFTs in $(**)$ are assumed to be periodic functions of period $N/2$ (where's rb-j when I need him?)  But, now note that $X_{\text{even}}[k]$ and $X_{\text{odd}}[k]$ are DFTs of length $N/2$, and so, if $N/2$ is an even number (hint: this means that $N$ is a multiple of $4$), we can use the idea behind $(**)$ to express $X_{\text{even}}[k]$ as
$$X_{\text{even}}[k] =  X_{\text{even, even}}[k] + W_{N/2}^kX_{\text{even, odd}}[k]$$
and $X_{\text{odd}}[k]$ as
$$X_{\text{odd}}[k] =  X_{\text{odd, even}}[k] + W_{N/2}^kX_{\text{odd, odd}}[k]$$
that is, both DFTs of length $N/2$ have been expressed in terms of 2 DFTs of length $N/4$, and so on.  If $N = 2^m$, we can proceed in this way to where we have $N/2$ DFTs of length $2$ that we compute using one multiplication and addition per DFT bin. So, $N$ multiplications and additions total to compute these $N/2$ DFTs. (Purists whining that a length $2$ DFT needs only one addition and one subtraction are kindly requested to go jump in the lake!). Next, we combine these $N/2$ DFTs of length $2$ into $N/4$ DFTs of length $4$ using one multiplication and one addition per DFT bin for a total of $N$ multiplications and additions to do the combination. Next, we combine these $N/4$ DFTs of length $4$ into $N/8$ DFTs of length $8$ using one multiplication and one addition per DFT bin for a total of $N$ multiplications and additions to do the combination..... Ultimately, we arrive at $(**)$.  This iterative procedure for computing the DFT of length $N$ is (one of the many versions of) the FFT algorithm for computing the Discrete Fourier transform. It has $m=\log_2N$ stages requiring $N$ multiplications and $N$ additions at each stage for a total of $Nm = N\log_2N$ multiplications and additions. In contrast, the use of Horner's rule would need $N(N-1)$ multiplications and $N^2$ additions. Hence the soubriquet FAST for this kind of algorithm.

Purists who have survived their jumps into the lake can now start their cacophony about radix-2 versus radix-4 versus radix-8 transforms and the exact number of multiplications and additions needed for implementing an FFT algorithm.


Estimate camera generalized velocity from consecutive frames

It's possible, but it takes some sophisticated processing to do it effectively in a general environment.

First of all, you can only detect motion relative to what you can see in the field of view. This requires good foreground/background separation to eliminate objects that might themselves be moving relative to the background (which is presumably what you want to measure the camera motion against).

Second, the background must have enough features (edges) in it to make all six degrees of freedom unambiguously observable.

It boils down to figuring out what the trajectories through the viewspace are of all the features in the scene, and then doing a cluster analysis of those trajectories in order to assign them to distinct objects. Finally, pick the cluster that most probably and most consistently represents the background. The inverse of this trajectory is the camera's trajectory.

Once you have the geometry of the situation squared away, you still need information about the size of objects in the field of view in order to scale the results correctly.


The most direct way to make these calculations is to design a set of filters to your specifications, then read the data out of the WAV files and run it through the filters. In Python/SciPy, you can use Python's wave module to read and write WAV files, SciPy's firwin to design a windowed sinc filter like SoX's sinc effect, and SciPy's lfilter to apply each filter. You will have to design each filter separately and call lfilter for each filter.

A simple way to approximate the energy in several bands is to first calculate the Discrete Fourier Transform (DFT) of consecutive frames, and sum up the energy in the bins corresponding to your bands of interest. The DFT approximates a signal with a sum of harmonically related sinusoids. The advantage of this technique is that it's easy to sum up several arbitrarily spaced bands. A big disadvantage is you have to apply a window to control spectral leakage. This makes the energy calculation an approximation. You can use NumPy's rfft function to calculate the Fourier Transforms.

If your bank of filters is related, for instance all the same bandwidth, or centered at harmonically or logarithmically related frequencies, there are some clever techniques that fall under the concept of filterbanks. This is more specialized, and I don't know of much pre-packaged code you could use.


Can Temperature Data be Predicted Using Adaptive Filter (Such As LMS) Algorithm?

Yes you can predict future temperatures, based on past temperatures, using adaptive filtering as well.

The optimal linear estimation of a WSS random process from its past values, which is known as linear prediction, is given by a Wiener filter structure where the desired response to be estimated is the current sample of the input (current temparature in your case) and the filter input is the $N$ past samples of the input, (assuming one step forward prediction of order $N$).

The LMS adaptive filtering algorithm simply approaches this optimal Wiener predictor coefficients for WSS signals and for non WSS signals tries to continue to be optimal by tracking it.

This prediction mechanism does not depend on the physical origin of the signals but on their statistical characterisation. As long as your temperature data posses reasonable degree of correlation within it, then the filter will do its best to predict it.


Scaling of the PWELCH function in MATLAB

I believe, there are two aspects of your question.

For the FFT calculation, the power density spectrum is the square of the FFT modulus multiplied by the sampling period (in your notation: T) and divided by the number of samples (in your notation:L):
T/L * abs(FFT(X))^2  or   T/L *abs(Y(f))^2

Here is now where I believe you have an error, you have to add the power of Y(f) and (Y(-f)) and not add first Y(f) and (Y(-f) and then perform the power calculation. By adding the amplitude first, in your notation,   P1(2:end-1) = 2*P1(2:end-1)  then after calculating the power, you have the double of the correct value.

A far as the Welch method for the calculation of the Power spectrum, it takes into account all the necessary normalization. You can chose between the power spectrum (option "power") or the power density spectrum (default option "psd")

Trying to use the option "power" for the welch calculation, it will gives as answer the amplitude of the sin squared divided by 2 (a sin as a power of 1/2): (A1^2)/2.

To calculate the power density spectrum, you have to divide the previous results by the applied frequency resolution, that is the sampling frequency by the number of samples, in your notation Fs/L. With the welch function, you have to use the option 'psd'.
Sorry for the long answer and hoping, this helps.


Why coherence function calculation needs averaging of different segments of the signal?

You can see how it falls out if you work through the simple case of N=2 (I've elided the 1/N for brevity)
\begin{equation}
\frac{(Y_{1}X_{1}^{*}+Y_{2}X_{2}^{*})(X_{1}Y_{1}^{*}+X_{2}Y_{2}^{*})}
{(X_{1}X_{1}^{*}+X_{2}X_{2}^{*})(Y_{1}Y_{1}^{*}+Y_{2}Y_{2}^{*})}
\end{equation}
It should be pretty obvious that there are now some cross products. Multiplying through:
\begin{equation}
\frac
{
X_{1}X_{1}^{*}Y_{1}Y_{1}^{*} + X_{1}^{*}X_{2}Y_{1}Y_{2}^{*} + X_{1}X_{2}^{*}Y_{2}Y_{1}^{*} + X_{2}X_{2}^{*}Y_{2}Y_{2}^{*}
}
{
X_{1}X_{1}^{*}Y_{1}Y_{1}^{*} + X_{1}X_{1}^{*}Y_{2}Y_{2}^{*} + X_{2}X_{2}^{*}Y_{1}Y_{1}^{*} + X_{2}X_{2}^{*}Y_{2}Y_{2}^{*}
}
\end{equation}
The second and third terms of the numerator are no longer the same as the second and third terms of the denominator.


If an image has large values at high frequency components then the data is changing rapidly on a short distance scale. e.g. a page of text

If the image has large low frequency components then the large scale features of the picture are more important. e.g. a single fairly simple object which occupies most of the image.

Thus low frequency correspond to slow varying information ( continuous surface)
high frequency correspond to quickly varying information ( eg. edges)


Received Morse Code signal-to-noise ratio?

Personally, I'd use whatever filtering scheme the actual receiver uses. If the narrow-band CW filter is exactly this, or even close, I'd use that. Then try to automatically (or manually) determine the regions where the signal is present. Measure the average power in the "off" regions, and that will be N. Measure the average power in the "on" regions, and that will be S+N. S = S+N - N, and then take 10*log10( S / N ). This will give you S/N ratio for the bandwidth of the narrow-band filter.


Matlab HIgh frequency sine wave generation

The function you need to evaluate is:

A*sin(2*pi*FREQUENCY*(1:NUM_SAMPLES)/SAMPLE_RATE)

Use a reasonable sample rate (at least twice the frequency) and a reasonable interval (remember that at a sample rate of 2Ghz, 1 second of signal will fill up the memory of your computer) ; and it'll work! From your code it seems that you are trying to generating 100s of signal at a sample rate of 1kHz, which doesn't make sense.

At the moment, in your code, the argument of the sin function is an integer multiple of $2\pi$ so in theory your code should be plotting a flat line of 0s. You are seeing something a bit different due to the finite resolution of floating point numbers.


Nyquist Rate (Sampling Frequency) for $ {f}^{2} \left( x, y \right) $

By the Convolution Theorem multiplication in Time / Spatial domain is equivalent of Convolution in the Frequency Domain.

The sampling rate (In its classic interpretation) is proportional to the support of a function in the frequency domain.

So if a function has a certain support in frequency, what would be its support after convolution with itself? Indeed it will be doubled in each dimension.


Newbie looking for an interesting introductory book to telecommunication and signals/data processing

"Telecommunications" is too broad a subject: you could spend your entire life learning about just a small part of the field. However, you mention data transmission specifically; I'll recommend these books:

"Digital Signal Transmission" by C. C. Bissell and D. A. Chapman, Cambridge. This is an undergraduate level book that has very lucid coverage of both the physical layer and larger telecom networks.

"Principles of Digital Communicatios", by B. Rimoldi, Cambridge. This is an advanced undergrad textbook that focuses on the math and theory behind digital communications. A great book if you like that kind of thing.

"Data and Computer Communications", by W. Stallings, Pearson. I'm not the biggest fan of Stallings, but his books are very accessible and cover all the bases. In fact, it might be a good idea to start with this book and then continue with more rigorous books like the first two I listed.


How to set sampling frequency of an acceleration data set

If you don't have a uniform time period (or nearly so) between samples, then you don't have a sample rate, and the FFT result may well be garbage.

If you don't know ahead of the sampling the  Nyquist frequency for the signal, then the sampling might be already contaminated by aliasing.

If the samples are close enough together relative to the highest frequency in the signal's spectrum, then you may be able to interpolate a new set of evenly spaced samples, which will then have a known sample rate (based on the distance of the new even spacing).


Interpretation of Convergence plot in terms of mean square error

0 dB corresponds to a MSE of 1. So, if your method gives MSE that approaches 0 dB as the number of iterations increases, then your method has an error floor; it cannot drive the MSE to any value smaller than 1 no matter how many iterations it takes (or how many data points you have) etc. A good method should drive the MSE to 0 (corresponding to $-\infty$ dB), that is, the MSE should decrease continually as the number of iterations increases, and not saturate at a value of 1 (0 dB). We can then argue as to whether the approach to 0 MSE ($-\infty$ dB) is exponential (MSE $= O(e^{-n})$) or polynomial (MSE $= O(n^{-k})$) including as special cases quadratic (MSE $=O(n^{-2})$) and
linear (MSE $= O(n^{-1})$), or sublinear
(e.g., MSE $= O(n^{-1/2})$ or MSE $= O((\ln n)^{-1})$ or similar)
where $n$ is the number of iterations (or number of data points).
Naturally, the faster the decrease, the better (and usually
more expensive) the method.

On the other hand, if the MSE of the best currently
known method is approaching
an error floor of 0 dB, there is a fantastic research opportunity (and
possibly a get-rich-quick opportunity too): devise a method that
drives the MSE to 0 instead of 1 and gain instant fame
(and if you patent the method and found a start-up company
to further develop it, possibly fortune RSN as well).


Low pass signal filtering using FFT or simple infinite impulse response filter

There are three potential problems with what you are doing:

Zeroing out the high frequencies causes ringing due to the Gibbs phenomenon.
You are essentially multiplying the Fourier transform with a square function, which is equivalent to circularly convolving the time-domain data with the inverse transform of the square function, which is a sinc function.  Unless you zero-pad the data before transforming it, the circular convolution will cause the data at the end of the FFT block to affect the data at the beginning, and vice versa.
If you were doing multiple FFT's to process a stream of data, unless you use a technique like Overlap-Add or Overlap-Save, there will be discontinuities at the boundaries between FFT blocks.  This problem is related to #2.


I and Q components and the difference between QPSK and 4QAM

Both QPSK and $4$-QAM constellations have signal points at $45, 135, 225$, and $315$ degrees (note typo in your question).  They arise from amplitude modulation (or, if
you prefer, phase modulation) of two carrier signals
(called the inphase and quadrature carriers) that are orthogonal (meaning that they
differ in phase by 90 degrees.  The canonical representation of a QPSK or $4$-QAM
signal during one symbol interval is
$$s(t) = (-1)^{b_I}\cos(2\pi f_c t) - (-1)^{b_Q}\sin(2\pi f_c t)$$
where $\cos(2\pi f_c t)$ and
$-\sin(2\pi f_c t)$ are the inphase and quadrature carrier signals at
frequency $f_c$ Hz
and $b_I, b_Q \in \{0,1\}$  are the two data bits (called the inphase
and quadrature data bits, naturally, since they are transmitted
on the inphase and quadrature carriers).  Notice that the
inphase carrier $\cos(2\pi f_c t)$ has amplitude $+1$ or
$-1$ according as the inphase data bit has value $0$ or $1$, and
similarly the quadrature carrier $-\sin(2\pi f_c t)$ has amplitude $+1$ or
$-1$ according as the quadrature data bit has value $0$ or $1$.
Some people regard this as an inversion of the normal scheme
of things, didactically asserting
that positive amplitudes must be associated with $1$ data bits and
negative amplitudes with $0$ bits. But if we look at it from the
phase modulation perspective, a $0$ bit means that the carrier
($\cos(2\pi f_c t)$ or $-\sin(2\pi f_c t)$ as the case may be) is
transmitted with no change in phase while a $1$ data bit creates
a change in phase (we will think of it as a phase delay)
of $180$ degrees or $\pi$ radians.  Indeed, another
way of expressing the QPSK/$4$-QAM signal is as
$$s(t) = \cos(2\pi f_c t - b_I\pi) - \sin(2\pi f_c t - b_Q\pi)$$
which makes the phase modulation viewpoint very clear. But,
regardless of which viewpoint we use, during a symbol interval,
the  QPSK/$4$-QAM signal is one of the following four signals:
$$\sqrt{2}\cos\left(2\pi f_c t + \frac{\pi}{4}\right),\\
\sqrt{2}\cos\left(2\pi f_c t + \frac{3\pi}{4}\right),\\
\sqrt{2}\cos\left(2\pi f_c t + \frac{5\pi}{4}\right),\\
\sqrt{2}\cos\left(2\pi f_c t + \frac{7\pi}{4}\right)$$
corresponding to $(b_I,b_Q) = (0,0), (1,0), (1,1), (0,1)$ respectively.

Note that the viewpoint taken here is of QPSK as consisting of
two BPSK signals on phase-orthogonal carriers. The demodulator
thus consists of two BPSK receivers (called the inphase branch
and quadrature branch, what else?).
An alternative view of QPSK as changing the phase
of a single carrier depending on a $4$-valued symbol is developed a little later.

The QPSK/$4$-QAM signal can also be expressed as
$$s(t) = \text{Re}\{B \exp(j2\pi f_c t)\}
= \text{Re}\{[(-1)^{b_I}+j(-1)^{b_Q}] \exp(j2\pi f_c t)\}$$
where $B$ is the complex-valued baseband symbol taking on values
in $\{\pm 1 \pm j\}$ and which, when plotted on the complex plane,
gives constellation points distant $\sqrt{2}$ from the
origin and at  $45, 135, 225$, and $315$ degrees
corresponding to data bits $(b_I,b_Q) = (0,0), (1,0), (1,1), (0,1)$
respectively. Note that complementary bit pairs lie diagonally across the circle from each other so that double bit errors are less likely than single bit errors. Note also that the bits naturally
occur around the circle in Gray code order; there is no need to massage a given data bit pair $(d_I,d_Q)$ (say $(0,1)$) from "natural representation" (where it means the integer $2 = d_I+2d_Q$: $d_I$ is the LSB and $d_Q$ the MSB here) to "Gray code representation" $(b_I,b_Q) = (1,1)$ of the integer $2$ as some implementations seem to insist
on doing. Indeed, such massaging leads to poorer BER performance since the decoded $(\hat{b}_I,\hat{b}_Q)$ must be ummassaged at the receiver into the decoded data bits $(\hat{d}_I,\hat{d}_Q)$ making the single channel  bit error $$(b_I,b_Q) = (1,1) \to (\hat{b}_I,\hat{b}_Q) = (1,0)$$
into the double data bit error
$${(d_I,d_Q) = (0,1) \to (b_I,b_Q) = (1,1)\to (\hat{b}_I,\hat{b}_Q) = (1,0)}\to (\hat{d}_I,\hat{d}_Q) = (1,0).$$

If we delay the four possible signals exhibited above by $45$ degrees or
$\pi/4$ radians (subtract $\pi/4$ radians from the argument of the cosinusoid), we get
$$\sqrt{2}\cos\left(2\pi f_c t + \frac{\pi}{4}\right)\Rightarrow \sqrt{2}\cos\left(2\pi f_c t + 0\frac{\pi}{2}\right) = \sqrt{2}\cos(2\pi f_c t),\\ \sqrt{2}\cos\left(2\pi f_c t + \frac{3\pi}{4}\right)\Rightarrow \sqrt{2}\cos\left(2\pi f_c t + 1\frac{\pi}{2}\right) = -\sqrt{2}\sin(2\pi f_c t),\\
\sqrt{2}\cos\left(2\pi f_c t + \frac{5\pi}{4}\right)
\Rightarrow \sqrt{2}\cos\left(2\pi f_c t + 2\frac{\pi}{2}\right) = -\sqrt{2}\cos(2\pi f_c t) \\
\sqrt{2}\cos\left(2\pi f_c t + \frac{7\pi}{4}\right)\Rightarrow \sqrt{2}\cos\left(2\pi f_c t + 3\frac{\pi}{2}\right) = \sqrt{2}\sin(2\pi f_c t),\\$$
which give the four constellation points at $0,90,180,270$ degrees referred
to by the OP.  This form gives us another way of viewing QPSK signaling:
a single carrier signal whose phase takes on four values depending on
the input symbol which takes on values $\{0,1,2,3\}$. We express this in tabular form.
$$
\begin{array}{|c|c|c|c|c|c|c|}
\hline
(b_I,b_Q) & \text{normal value} ~k & \text{Gray code value} ~\ell & \text{signal as above} &\text{phase-modulated signal}\\
\hline
(0,0) & 0 & 0 & \sqrt{2}\cos(2\pi f_c t)
& \sqrt{2}\cos\left(2\pi f_c t - 0\frac{\pi}{2}\right)\\
(0,1) & 1 & 1 & \sqrt{2}\sin(2\pi f_c t)
& \sqrt{2}\cos\left(2\pi f_c t - 1\frac{\pi}{2}\right)\\
(1,1) & 3 & 2 & -\sqrt{2}\cos(2\pi f_c t)
& \sqrt{2}\cos\left(2\pi f_c t - 2\frac{\pi}{2}\right)\\
(1,0) & 2 & 3 & -\sqrt{2}\sin(2\pi f_c t)
& \sqrt{2}\cos\left(2\pi f_c t - 3\frac{\pi}{2}\right)\\
\hline
\end{array}
$$
That is, we can regard the QPSK modulator as having input
$(b_I,b_Q)$ that it regards as the Gray code representation
of the integer $\ell \in \{0,1,2,3\}$ and produces the
output
$$\sqrt{2}\cos\left(2\pi f_c t - \ell\frac{\pi}{2}\right).$$
In other words, the phase of carrier $\sqrt{2}\cos(2\pi f_c t)$ is
modulated (changed from $0$ to $\ell\frac{\pi}{2}$) in
response to the input $\ell$.

So how does this work in real life or MATLAB, whichever comes first?
If we define a QPSK signal as having value $\sqrt{2}\cos\left(2\pi f_c t - \ell\frac{\pi}{2}\right)$ where the value of $\ell$ is typed in as 0 or 1 or 2
or 3, we will get the QPSK signal described  above, but the
demodulator will produce the bit pair $(b_I, b_Q)$ and we must
remember that the output is $\ell$ in Gray code interpretation,
that is, the demodulator output will be $(1,1)$ if $\ell$ happened
to have value $2$, and interpreting output $(1,1)$ as $3$ is a decoding
error that is not generally discussed in textbooks!


stationary vs. undecimated wavelet transform

From what I understand un-decimated wavelet transform (UDWT) involves the elimination of the decimation step to enable DWT to have translational invariance property. The stationary wavelet transform is one such algorithm which is also known as algorithme à trous. There are a number of other algorithms of UDWT a quick google scholar search can you show that but they all use similar ideas to the stationary wavelet transform.


MATLAB's $\tt dsp.ParametricEQFilter$ bandwidth question

Apparently, after some research, I have found that the bandwidth "bandedges" is a bit of a tricky subject: (see Bandwidth Re-visited).
Some references suggest indeed that the bandwidth is measured at the gain which is equal to half of the peak gain. Some others say that you should measure the bandwidth at $ +3\ \textrm{dB}$ for boost gain, $ -3\ \textrm{dB}$ for cut. (which eventually runs into the issue of where it should be measured if the gain is less than $ 3\ \textrm{dB}$.

Hilmar is right about the geometric mean of the corner frequencies. However, the cutoff points (i.e. the bandedges) do not exist at half the peak gain in MATLAB.   They exist at $G^{2}_B = \frac{1+G^2}{2}$, where $G$ is the parametric equalizer linear gain, and $G_B$ is the linear bandwidth gain as described here.


Why does the impulse response determine the transfer function of a system?

Because every signal can be decomposed into a linear series of scaled impulses that are shifted in time. Thus applying a linear time-invariant system on each impulse in the series and then summing the results again will give the same result as applying it on the whole signal. Therefore, one only needs to know how the system responds to an impulse to be able to calculate its response to an arbitrary signal and all properties of the system.
This DSP guide chapter has some nice illustrations on the subject.
This DSP guide chapter (only in PDF) shows how from the impulse response as written in the form of a difference equation the transfer function can be derived using the z-transform. For FIR filters the $b_i$ coefficients will be zero in equation 33-2.


How does Siri recognize me saying "Hey Siri"?

"Ok Google" is described in many publications by Google

Automatic Gain Control and Multi-style Training for Robust Small-Footprint Keyword Spotting with Deep Neural Networks

Convolutional Neural Networks for Small-Footprint Keyword Spotting

It is based on DNN specifically trained for keyphrase and runs really fast. It does not consume a lot of power even on mobile.

Alexa spotting is implemented by the same guys and available as Snowboy

UPDATE: Apple describes their implementation here.


low pass filter and FFT for beginners with Python

The fact that the result is complex is to be expected. I want to point out a couple things:

You are applying a brick-wall frequency-domain filter to the data, attempting to zero out all FFT outputs that correspond to a frequency greater than 0.005 Hz, then inverse-transforming to get a time-domain signal again. In order for the result to be real, then the input to the inverse FFT must be conjugate symmetric. This means that for a length-$N$ FFT,

$$
X[k] = X^*[N-k], k = 1, 2, \ldots , \frac{N}{2} - 1\;\;\;\;\;\;\;(N\;\; even)
$$

$$
X[k] = X^*[N-k], k = 1, 2, \ldots , \lfloor \frac{N}{2} \rfloor\;\;\;\;\;\;\;\;(N\;\; odd)
$$

Note that for $N$ even, $X[0]$ and $X[\frac{N}{2}]$ are not equal in general, but they are both real. For odd $N$, $X[0]$ must be real.

I see that you attempted to do something like this in your code above, but it is not quite correct. If you enforce the above condition on the signal that you pass to the inverse FFT, then you should get a real signal out.

My second point is more of a philosophical one: what you're doing will work, in that it will suppress the frequency-domain content that you don't want. However, this is not typically the way a lowpass filter would be implemented in practice. As I mentioned before, what you're doing is essentially applying a filter that has a brick-wall (i.e. perfectly rectangular) magnitude response. The impulse response of such a filter has a $sinc(x)$ shape. Since multiplication in the frequency domain is equivalent to (in the case of using the DFT, circular) convolution in the time domain, this operation is equivalent to convolving the time domain signal with a $sinc$ function.

Why is this a problem? Recall what the $sinc$ function looks like in the time domain (below image shamelessly borrowed from Wikipedia):

The $sinc$ function has very broad support in the time domain; it decays very slowly as you move in time away from its main lobe. For many applications, this is not a desirable property; when you convolve a signal with a $sinc$, the effects of the slowly-decaying sidelobes will often be apparent in the time-domain form of the filtered output signal. This sort of effect is often referred to as ringing. If you know what you're doing, there are some instances where this type of filtering might be appropriate, but in the general case, it's not what you want.

There are more practical means of applying lowpass filters, both in the time and frequency domains. Finite impulse response and infinite impulse response filters can be applied directly using their difference equation representation. Or, if your filter has a sufficiently-long impulse response, you can often obtain performance benefits using fast convolution techniques based on the FFT (applying the filter by multiplying in the frequency domain instead of convolution in the time domain), like the overlap-save and overlap-add methods.


Fourier series of a square pulse (simple simplification)

You could rewrite $x(t)$ as a rectangular function before doing the Fourier transform. You could deduce the width and center of the rect signal directly from the given expression but it would be a easier to just draw the signal and deduce the values from the graph.

Your rect signal will be in the following form:

$x(t) = rect(\frac{t-t_o}{W})$

where:

$t_0$ - center of the rectangular signal

W - width(duration) of the pulse

From this point you would use a Fourier transform table and some Fourier transform theorems (involving scaling and delay) to get the $sinc$ version of the $X(\omega)$.

Or, if you prefer, you could just modify a little bit the expression you already have to get the expression where you could use Euler's formula as Dilip Sarwate commented.

$X(\omega)=1 \cdot \frac{1}{j\omega}(e^{-j\omega} - e^{-j2\omega})= e^{-j1.5\omega}e^{+j1.5\omega} \frac{1}{j\omega}(e^{-j\omega} - e^{-j2\omega}) = e^{-j1.5\omega}\frac{1}{j\omega}(e^{+j0.5\omega} - e^{-j0.5\omega})$

From this point it should be easy to get the result.


How we can encode/decode sparse signals?

Since you have tagged your question compressive-sensing, I would guess that you are interested in a compressive sensing solution to this. In principle, this could work by taking a number (less than 64) of random linear combinations (measurements) of your length 64 vector and transmitting these. At the receiver, the fact that this vector has a sparse representation allows you to reconstruct it again from these measurements. However, compressive sensing does not work that well on low-dimensional vectors (64 is relatively small) and so I do not expect this scheme to work very well.

A more natural solution, since you mention that you know the sparse length 128 representation of the vector, is to use run-length encoding (Wikipedia explanation) of the sparse representation and transmit that. Since the sparse vector is almost all zeros, this should be very efficient.


Better way to mute fundamental baseband frequency?

Use a notch filter such as the one shown in the figure below, where $\omega_n =2\pi 60/f_s$, where $f_s$ is your sampling rate, and $\alpha$ is chosen based on the bandwidth of the notch and how long you can allow for settling in the time domain; the tighter the filter bandwidth the longer it will take to settle; you can use a first order approximation of 10%/90% rise time = $\frac{0.35}{BW}$ where the BW is in Hz (and the amount of precision you can allow in the filter is an important factor---insufficient precision and $\alpha$ to close to 1 can lead to instability).

See this post for details: Transfer function of second order notch filter

Note the diagram below to facilitate implementation, mapping the coefficients to the transfer function as commonly referenced in Matlab and Octave. Do note that the signal flow below is going from right to left, which is opposite the diagram above (this was because my slide was showing how the implementation came from transposing the Direct Form II implementation). To be clear, the resulting coefficients are:

$b_0 = 1$

$b_1 = -2cos(\omega_n)$

$b_2 = 1$

$a_1 = -2\alpha cos(\omega_n)$

$a_2 = -\alpha^2$

A normalizing gain of $(1+\alpha)/2$ is added to the whole filter. Since $\alpha$ is typically close to 1, this can reasonably be omitted in most cases.


How do I plot the spectrogram of a wav file?

This isn't a great question, but I'll try to answer:

  All the files that I have, have two channels.

A spectrogram can only show one channel at a time.  Either make one spectrogram for each channel, or average your channels together into one.  Here's an example of the syntax I found with a quick google search:

[y, fs]=wavread(fileName); % Read wave file
left=y(:,1); % Left channel
right=y(:,2); % Right channel

then you can do spectrogram(left)

  I know that the Spectrogram is supposed to give you an idea of the fundamental frequencies present in your audio file.

Not quite.  The term "fundamental frequency" means something different.

  What does the Spectrogram exactly do?

It shows the frequency spectrum of the signal and how it changes over time.  If the beginning of the signal is white noise (flat spectrum), and the end of the signal is a tone (spikey spectrum), the spectrogram will show how it changed from one spectrum to the other over time.  It does this by dividing the signal up into small chunks and calculating the spectrum of each chunk.

  What do the parameters mean? This is the structure of the Spectrogram function call : S=spectrogram(x,window,noverlap,nfft,fs). What do the window and noverlap parameters mean?

Did you read the documentation?

window is a window function that is applied to each chunk of the file before the frequency transform.  If you didn't use a window function, there would be abrupt changes at the beginning and end of each chunk, which produce wide frequency components that you don't want.  The window function is like a "fade-in, fade-out" envelope, so that you only see the spectrum of the middle of the chunk, and not the artifacts caused by chopping the chunk abruptly at the endpoints.

noverlap is the number of samples that each chunk overlaps the next.  The windowing function attenuates the signal to near 0 at the boundaries of each chunk, so if you didn't overlap the chunks at all, you would be emphasizing parts of the signal and ignoring other parts.  The overlap lets you measure the same part of the signal twice, once in the middle of the chunk where it's important, and once at the ends where it's mostly ignored.

  How do I set these parameters?

Basically you try different ones and see which is best at showing you what you want to see.  Different window functions are designed to do different things, so read about them and what they are meant to do.

  Are they input file dependent?

Only if you're looking for different types of things in different files.


By projecting a vector x using PCA (on the PCs), you maximize the variance in the reduced space. Initially, the space is not optimal in terms of maximizing the variance.

So:

  PCA projects vector  to a space of  dimensions where the
  difference between the initial vector and the projection has maximum
  energy.

(initially the is no maximum variance but after projection the space is optimal in terms of maximizing the variance in the PC defined space)

PS:

Minimizing J1 (the error) corresponds to maximizing the quadratic form e'Se. After demeaning this is e'X'Xe. Set now w = Xe (the projection of the data X). The variance in the projected space is w'w that is equal to the e'X'Xe. Thus, minimizing the error means maximizing varinace


Range is the minimum to the maximum value of something. Since bandwidth is the difference between the upper and lower frequencies in a continuous band of frequencies, then the range of bandwidth of 15 to 20 KHz would then mean that the minimum difference (bandwidth) would be 15 KHz and the maximum bandwidth would be 20 KHz. If the saying was instead that the range of occupied frequencies was 15 to 20 KHz then it would then mean the bandwidth was 5 KHz.


You can analyse the input online, and get some rough estimates on its power, energy, pdf, autocorrelation matrix, etc. and also its eigenvalues of (time-varying) autocorrelation matrix. Then you may estimate the associated upper bound for the step-size of LMS.
Nevertheless, finding the optimal step-size in applications within unknown environments is not easy. So many of the LMS variations are based on different approaches for efficiently, and reliably finding the best step-sizes.


What frequencies are present in the Fourier transform of the Dirac impulse?

A Dirac impulse $x(t)=\delta(t-d)$ has the continuous-time Fourier transform $X(\Omega)$ of
$$\mathcal{F}\{\delta(t-d) \} = 1 e^{-j\Omega d} $$

whose magnitude is $$|X(\Omega)| = 1 ~~~, \text{ for all } \Omega $$ and a phase of $$<X(\Omega) = -\Omega \cdot d $$

Note that it's incomplete to think of the real or imaginary parts of the Fourier transform alone. Rather the magnitude and phase point of views are more reflective of the nature of the result.

So in this case the magnitude is $1$ and hence it's said to contain all frequencies of magnitude $1$.
Note that these are differential amplitude components of continuum frequency range as opposed to a finite amplitude of discrete set of frequency components, aka line components.


Using Matlab's filter() function for solving difference equation with initial conditions

The problem is the definition of the initial condition of the function filter(). Note that according to the mathworks reference page, filter is implemented using a direct-form II transposed structure. If you look at the corresponding diagram, you see that the filter states (called $w_k(m)$ in the figure) are defined after multiplication with the respective filter coefficients.

So in your example if you want to initialize filter with some past output value $y[-1]$, your initial state is given by $y[-1]\cdot 0.5$ because that's the value of the input of the delay element.

y = filter(num_coef, den_coef, x, 0.5)

should give you the expected result.

Note that for the general case there's the function filtic, which computes the necessary initial conditions given past output samples (and, if necessary, past input samples). For your example you would use the command

zi = filtic(num_coef, den_coef, 1)

which results in zi = 0.5, as explained above.


Meaning of negative frequencies in Baseband non sinusoidal, non periodic signal

It means the same thing.
The negative frequencies are there to cancel out the complex portion of the signal in the time domain, so that the time domain signal is constrained to the real axis.
The Fourier Transform decomposes a signal, $x(t)$ onto a set of complex basis functions: $e^{-j\omega t}$ for $\omega \in (-\infty, \infty)$.
If $x(t)$ is real, then the decomposition, i.e. the Fourier Transform $X(\omega)$, must be such that the complex portions cancel out, when the inverse transform is performed to recompose the signal from it's constituent parts.
This is a conceptual explanation of why the Fourier Transform of a real signal is symmetric about $0$ Hz.


Since you have defined your signal $x[n]$ to be a set of non-negative numbers with a single global maxima, the algorithm to find the maximum is quite simple.

compute the first (causal , backward) difference
$$y[n] = x[n] - x[n-1]$$

Then look for the sign change of $y[n]$; i.e., at the extreme point you will have the sign of $y[n]$ changing from positive to negative,  assuming you had a single global maxima.

$$z[n] = \text{sign_change_detector} \{ y[n] \} $$

Don't forget to check the end points as well. As they are also the candidates to be maximum in addition to those critical point(s).


Laplacian

It is a 2nd derivative of the recorded voltage for each electrode. Suppose the recordings at position $(i,j)$ is $V(i,j)$, then Laplacian operation  at $(i,j)$ is

$(V(i,j)-2V(i-1,j)+V(i-2,j))^2 + (V(i,j)-2V(i,j-1)+V(i,j-2))^2$

And it serves as a high-pass filter that enhances localized activity while suppresses the diffusion activity. The implementation of the Laplacian in EEG filtering on the voltage at each electrode is to subtract the weighted voltages from the surrounding electrodes from the voltage recording at current electrode, where the weight is electrode distance dependent.

Large

The difference between large and small is determined by how you select the electrodes. For the small Laplacian, you select the nearest-neighbor electrodes, while for the large Laplacian, they are the next-nearest-neighbor electrodes. See the figure below. If the signal you are interested in is less highly localized and/or varies in exact location over
time, the large Laplacian is expected.

Read the paper Spatial filter selection for EEG-based communication by Dennis J. McFarland et all for more details.


Does a signal and its autocorrelation ever have the same Fourier spectra?

(a). Yes, a signal and its autocorrelation have the same frequencies. If $X(f)$ is the Fourier transform of the signal and $A(f)$ the Fourier transform of its autocorrelation, then $A(f) = |X(f)|^2$. So, in general, a signal and its autocorrelation don’t have the same Fourier transform, but there is at least one exceptional case:  $x(t) = W \operatorname{sinc}(Wt)$ whose Fourier transform is $X(f) = \operatorname{rect}\left(\frac fW\right)$ and so $A(f) = |X(f)|^2 = X(f)$. Those who wish to generalize this notion to signals whose Fourier spectra are the sum of non-overlapping $\operatorname{rect}$ functions should write out the details for themselves.

All of the above is for finite-energy signals.

If one is thinking of periodic finite-power signals with Fourier series $x(t) = \sum_{n=-\infty}^\infty c_n \exp(j2\pi nf_0t)$, then the periodic autocorrelation function has Fourier series $\sum_{n=-\infty}^\infty |c_n|^2 \exp(j2\pi nf_0t)$. Fourier series are represented by impulse trains in the frequency domain:
\begin{align}
X(f) &=  \sum_{n=-\infty}^\infty c_n \delta (f-nf_0),\\
A(f) &= \sum_{n=-\infty}^\infty |c_n|^2 \delta (f-nf_0)
\end{align}
and if we wish to continue to use the formula $A(f) = |X(f)|^2$ with the spectra of these periodic signals, then we must assume that
$$c_n \delta(f- nf_0)c_m^* \delta^*(f- mf_0) = \begin{cases}
|c_n|^2 \delta(f- nf_0), & \text{if}~m=n,\\
0, &\text{if}~m\neq n.
\end{cases}$$
Note that we are assuming that $|\delta(f-nf_0)|^2 = \delta(f-nf_0)$ and $\delta(f- nf_0) \delta^*(f- mf_0)=0$ for $m\neq n$, but if we are willing to swallow this bald and unconvincing narrative, then everything is hunky-dory.  Once again, we see that the signal and its autocorrelation have the same frequencies but that their Fourier transforms are generally not equal unless it so happens that $|c_n|^2 = c_n$ for all $n$. Note that this means that for each $n$, $c_n$ must equal $1$ or $0$. Thus, the only real-valued periodic signals whose Fourier transforms equal the Fourier transform of their autocorrelation function are of the form
$a_0 + \sum_{n=1}^\infty a_n 2\cos(2\pi f_0 t),~~ a_i \in \{0,1\}, i = 0, 1, 2, \ldots$


Remove brighter regions in an image on frequency domain

It is important to understand what the frequency domain means. The bright regions in an image correspond to pixels with a certain value, whether it be color or gray scale. Frequency domain has the interpretation of looking at how fast the values change, not their values. For example, think about an image of a dark blue car against a white background. The pixels along the edge may change from dark blue to white just by moving a few pixels, this is high frequency because the change in pixel value is high over a just a few pixels. This is different from a picture where you have dark blue which slowly fades to white over many pixels. Filtering in the frequency domain on an image would be affecting how the pixels change in value. For example, using a low pass filter on the dark blue car picture I mentioned, would block the high frequency components. This means the sharp changes in pixel value, like on the edge between the dark blue car and white background, would be smoothed (blurred) out.

Is it possible to remove white pixels from an image with frequency domain filter? No, and if you need to do that it would be better to just scan the image for pixels with a white value or value above some threshold and remove those pixels.


Your fft vector Y has the same length as your input signal y, because you just specified nf without letting fft() know what your desired FFT-length is. This is why your peak does not appear at 440 Hz. Your vector Y does not correspond to the frequencies in f. It's just a matter of correct scaling of the x-axis.

EDIT: I do not know scilab, so I don't know how to pass the desired FFT-length to fft(). If you can't do that, you just need to make your time domain signal have the desired FFT-length.


The optimal decision regions are the Voronoi Regions. I dont know, if this is what you are looking after.

import numpy as np

points = np.array([(1,1), (1,-1), (-1,1), (-1,-1), (3,3), (3,0), (3,-3), (0,-3), (-3,-3), (-3,0), (-3,3), (0,3), (5,0), (0,5), (-5,0), (0,-5)])
from scipy.spatial import Voronoi, voronoi_plot_2d
vor = Voronoi(points)
voronoi_plot_2d(vor)
plt.axis('equal'); plt.xlim((-7,7)); plt.ylim((-7,7));

I actually, consider analytically calculating the error probability of this constellation a very hard task, since the regions are not very regular.


Echo generation of FIR system using Matlab

Generally, an echo of a signal is not just its delayed and weakened copy, but is characterized by an impulse response of the echo path. The basic model of an echo is a linear one.
In order to simulate a linear echo $E$ one will implement the following convolution:
$E = S * H, $
where $S$ is your sampled signal and $H$ is an echo path impulse response, both with the same sampling frequency.
The particular form of the impulse response $H$ depends on the echo path. For example, $H$ may include $N$ zero values at the beginning to simulate an echo delay $N$.
Note, that in case of an acoustic echo, the linear impulse response $H$ may be not enough since an echo path is usually non-linear.


Frequency-domain deconvolution: "Direct" filtering vs "Wiener" filtering

The direct approach is noise sensitive and the second one (which is also known as $H_1$ estimator) is somehow noise resistant.
$H_1$ estimator assumes that there is no noise at the input signal, $B(\omega)$ in this case, and the noise exists only at output $D(\omega)$ and is uncorrelated to input $B(\omega)$. Therefore, the observed output signal is
$$
D'(\omega) = D(\omega)+N(\omega)
$$
where $N(\omega)$ is the uncorrelated noise. Under this assumption, the $H_1$ estimator gives a better result
$$
F'(\omega) = \frac{D'(\omega)B^*(\omega)}{B(\omega)B^*(\omega)} = \frac{D(\omega)B^*(\omega)+N(\omega)B^*(\omega)}{B(\omega)B^*(\omega)} = F(\omega) + \frac{N(\omega)B^*(\omega)}{B(\omega)B^*(\omega)}
$$
Since the noise is uncorrelated to the input signal, the second term can be eliminated by time averaging.

Here is some code and you can try some other SNR values. In practice for real data, you cannot obtain the clean signal without noise. All signal you measure is kind of noisy.
clear

% Input signal
x = mls;

% Comb filter (unknown system)
b = [1, zeros(1, 50), 1];
y = filter(b, 1, x);
fvtool(b, 1)

% Add white noise to the output
SNR = 10; % try some lower value and see what's happening
yn = awgn(y, SNR);

% H1 estimator
pxy = cpsd(x, yn);
pxx = cpsd(x, x);
H1 = pxy./pxx;
figure; subplot(211); plot(mag2db(abs(H1)))

% direct division
X = fft(x);
Y = fft(yn);
X = X(1:ceil(length(X)/2));
Y = Y(1:ceil(length(Y)/2));
H0 = Y./X;
subplot(212); plot(mag2db(abs(H0)))


Physical meaning of average values of random signals

Dilip: I'm confident that what you said is correct but the weirdo sample path that you refer to ( the one with probability zero ) can be used to estimate the autocorrelation at different lags, namely $\hat{R}_{xx}(\tau)$. This is true even though it's only one realization and this is the case only because of the ergodicity assumption. So, to a person coming in new, my opinion is that your answer doesn't give enough significance to the one realization. It's a weirdo but a useful weirdo for sure  !!!!!


Using the system function $H(z)$ to the solve the output of a DT LTI system

Let's first briefly recall the meaning of $H(z)$.
Let $y[n]$ be the output of the LTI system. Given an input $x[n]$, the output can be written
$$y[n] = h[n] *x[n] = \sum_{k=-\infty}^\infty h[k]x[n-k]$$
where $h[n]$ is the impulse response of the LTI system en $*$ represents the convolution. Let's now assume that $x[n]$, the input of the LTI system, is of the form $z^n$. In this case,
$$y[n] = \sum_{k=-\infty}^\infty h[k]z^{n-k} = z^n\left(\sum_{k=-\infty}^\infty h[k]z^{-k}\right) = z^nH(z)$$
where $H(z)$ denotes the transfer function of the LTI system (the z-transform of the impulse response). It can be seen that that $z^n = e^{j\Omega}$ is an eigenfunction of the system, that is, it is only affected by a complex factor (i.e. $H(z)$, the so-called eigenvalue) at the output of your system.

Now, your input $f[n]$ can be written as
$$f[n] = \mathcal{R}\{Ae^{j(\frac{\pi}{2}n + \frac{2\pi}{10})}\} = \mathcal{R}\{Ae^{j\frac{2\pi}{10}}e^{j\frac{\pi}{2}n}\}$$
and $Ae^{j\frac{2\pi}{10}}e^{j\frac{\pi}{2}n}$ matches the form of an eigenfunction (with a complex factor $Ae^{j\frac{2\pi}{10}}$). The output of your system to $Ae^{j\frac{2\pi}{10}}e^{j\frac{\pi}{2}n}$ is then given, using linearity of your system and the polar notation of $H$, by
$$y[n] = Ae^{j\frac{2\pi}{10}}H(e^{j\frac{\pi}{2}})e^{j\frac{\pi}{2}n} = Ae^{j\frac{2\pi}{10}}|H(e^{j\frac{\pi}{2}})|e^{\arg\{H(e^{j\frac{\pi}{2}})\}}e^{j\frac{\pi}{2}n}.$$
Finally, again by linearity, the output of your system to $f[n]$ is given by the real part of that.