So the 100Hz sine wave is sampled at 2000Hz resulting in a fft output mirror about 1000Hz, the nyquist frequency.

no. The FFT doesn't give you that. It's the spectrum as visualized here. It's an interpretation of the FFT (which is just an implementation of the DFT). For the DFT, it's just a single periodic function with a period of 5. That simple. No repetitions.

The repetitions in the visualization stem from the fact that mathematically, discrete spectra are periodic, and that the spectrum graphs were chosen to have frequency axis longer than a single Nyquist bandwidth. Again, the DFT doesn't care about what it means – it will just "show" that there is a periodicy of frequency $\frac 15$.

  It is then downsampled, effectively sampling this signal at 500Hz (factor of 4)/ or is it more a appropriate to say the signal now has characteristics like the 100Hz wave was sampled at 500Hz? This is achieved by discarding samples.

You're starting to get into the right mindset. Yes, this is like a 100 Hz signal sampled at 500 Hz – count the samples! There's a repetition every five samples!

Generally, I recommend just "dropping" knowledge of the original frequency at the ADC. Really, it makes things easier. You're sampling a 100 Hz tone with 2000 Hz, so your period is 20 samples, and the relative frequency is $\frac 1{20}$. Then you decimate by 4 — leading to a relative frequency of $4\frac1{20}=\frac15$.

The DSP algorithms don't have any of that "original frequency" knowledge; as said, the DFT really doesn't care whether it's a 100 Hz signal sampled at 2 kHz, or it's a 4 MHz signal sampled at 80 MHz – all that counts is how long a period is, in samples, because all the DSP only sees your signal as sequence of numbers. One after the other – no matter how fast the orignal samples came in.

  This wave is then upsampled by four, however I do not see what this does to benefit the signal?

Wrong question to ask. It's done to illustrative interpolation.

  I understand afer this it is low-pass filtered to remove the new alias images due to downsampling, but surely you you could low-pass filter after the downsampling and still end up with the same frequency output, i.e. 100Hz?

Nononononono!

You didn't learn your lesson; yes, your signal after decimation could be filtered, but that would have no effect, because it still fulfills the bandwidth criteria set by Nyquist. That's why there's no anti-aliasing filter here – the signal is already sufficiently band-limited as is.

However, after interpolation, you will need an anti-imaging filter.

Because the inserting of zeros did exactly that, adding images where there should be none.

  Or does it allow a 100Hz output with increased data points, which is better?

Again, not the question. This is an illustration, and the fact that you wanted to omit the anti-image filter shows that it's a good illustration: compare $|X_i|$ and $|X_u|$. You only want to see what's left in $X_i$, not all the spectrum that is in $X_u$!


short time energy matlab code for audio feature extraction

Short time energy (STE) is simple feature and so is the code to compute it.

ste = sum(buffer(signal.^2, winLen));

where  and winLen is duration of rectangular time window (in samples).

If you want to use diffrent shape of window, then you have to multiply each column returned by buffer by this window, before you sum the columns.


The most practical attempt that I am aware of is by Won and Berger (2005). They simultaneously recorded vocalizations at the mouth with a microphone and on the skull with a homemade vibrometer. They then estimated the relevant transfer functions with linear predictive coding and cepstral smoothing.


Idea for Noise Level Estimation / Automatic Thresholding in the Presence of Peaks

I'd do some small adjustments to your idea (You really nailed them).

Assumptions

The Signal Model - Signal + Additive White Gaussian Noise (AWGN)
Probably we could generalize it more but this is beyond the scope of this question.
The DFT of the signal contains Peaks with relatively small roll off
This is important as we're almost saying the Signal is a combination of sparse number of Harmonic signals. As the lobes width will hide the noise data. So the ratio between the observation window, signal lobes width and sampling rate means that the DFT of the data creates sharp peaks.

Estimation

Under the above assumptions we know that:

Noise is AWGN in the DFT (As Gaussian Vector stays Gaussian under linear transformations). It also suggests that the Mean Value and the Median are similar (Realization of Symmetric Distribution).
Without the Signal Data the Mean and the Median value of the beans should be similar.
Number of bins contaminated with signal (Funny to write this, ha?) is sparse.

So we can sort the bins, but sorting and evaluating the mean on sub set will create bias unless we make sure the sub set is symmetric around the Mean / Median.
So, I'd sort the data and keep bins in the 25-75 Percentile range only. I'd use Mean on those.


In general, Maths to audio signal processing is a very natural progression. I seem to recall my Master's and PhD supervisors both did their undergraduate degrees in Maths.

To your questions:

This is the hardest question to answer. Depends what's in your "portfolio of basic competence". In addition to the understanding of signal processing theory that you're developing you need to show you're able to implement these algorithms probably either in hardware or software. Given little EE grounding, I'd suggest software. Look at developing your software engineering skills. They will serve you well in either academia or industry.
MSc will definitely help but is not absolutely necessary. The types of algorithms you are describing could conceivably be the subject of an MSc level project. I work at a company that makes hardware for audio signal processing. The statistic they gave me is that around 60% of staff have a Master's or higher degree.
Yes, combining signal processing and machine learning is a very wise idea if you want to work towards the research end of the field. Several problems in audio signal processing - e.g. audio source separation and speech recognition- are increasingly finding their state-of-the-art solutions come from machine learning.


Why is random noise amplitude inversely proporional to the frequency?

This looks like a pink noise model, where the power spectral density of the noise is proportional to $\frac{1}{f}$. This model is appropriate for some real-world phenomena; I can't speak to whether it's the right choice for this particular application.

Your interpretation of the meaning of this sort of spectrum is correct, although the noise is still considered random. Instead, we consider consecutive noise samples to be correlated with one another (in contrast to white noise, where all samples are jointly uncorrelated). That doesn't mean that the signal isn't random.


Is it applicable to implement a cascaded moving average filter with variable window on real-time?

yep, that's totally doable.
However, since you want to be able to adapt at any point, the simple "add newest, subtract oldest sample from sum" approach to moving averages might no longer work (you'd need to be able to subtract or add more or fewer samples on the fly, which might require you to recalculate the whole sum, within one sample). At that point, the MA loses its computational advantage.

If so, how can I switch between different window sizes without causing spikes in filtered data or abrupt changes in the rate of the smoothed data?

Interpolate between the old and the new filter impulse response, and cycle through the interpolated filters slowly enough for you to avoid what you define to be a spike. You can't mathematically have "no abrupt change" but "abruptly change". That's the point where your Moving Average really needs to be treated as any other direct form FIR filter.
So: why a moving average? Unless you're really looking for a filter with a sinc spectrum, that's rarely the most adequate smoothing filter. And through your interpolation requirements, you lose most of the computational advantage of a simple moving average implementation.


Find original Hz of converted audio-file?

The short answer is no.

The longer answer is that if there is some known signal that you know the frequency you could but this is knowledge you bring to the problem.

You can look around 60Hz or 50Hz depending what the utility power frequency uses and see it there is some suggestion of a steady line in the data in that frequency range.

Musicians, after thousands of hours of playing an instrument will know what it should sound like.  While somewhat subjective, musical aesthetics are more than adequate for aesthetic purposes.

Summary,
If all you have is a file, not a chance.

If you have some clues and hunches and familiarity with the data there is a reasonable possibility.


State of the art algorithms for estimating noise variance

If you hope that $f[n]$ can be sparser in some linear representation $L$, then robust statistics can be used to estimate variance bounds on coefficients of $L(X[n])$. One common method is orthogonal wavelet denoising, where a common standard deviation estimate is:
$$ \mathrm{median} \frac{|c_i|}{0.67449}
$$
where $c_i$ denotes the highest scale subband. The magic $0.67449$ scale factor relates the standard deviation to the robust Median absolute deviation through the error function erf ("MAD equals the half-normal distribution median").

How robust is that? I tested it on a classical "crane" image, with several noise realizations, and here is the $\sigma$ versus  $\hat{\sigma}$ estimates with dispersion.

Not to bad: somehow linear, albeit with an offset at zero, and an increasing  dispersion. This can suffice for some uses. Meanwhile, for a better interval of estimates, I usually use several representations (wavelets, filter banks) to capture unknown sparsity optimality, with a certain amount of redundancy, and combine them into a range of estimates. Since redundancy kills independence, a calibration is required to get thre proper factor, that easily done with transforming white noise with the set of transformations you have. This is what we use with redundant Optimization of Synthesis Oversampled Complex Filter Banks or for denoising with Noise Covariance Properties in Dual-Tree Wavelet Decompositions.

Although the estimate is excellent with simulated Gaussian noise, this is of course less efficient in practice, as noise estimated from one realization is rarely Gaussian (or something easy to handle).

Additional references for images:
- Noise Estimation from a Single Image, with Matlab code

  In order to work well, many computer vision algorithms require that
  their parameters be adjusted according to the image noise level,
  making it an important quantity to estimate. We show how to estimate
  an upper bound on the noise level from a single image based on a
  piecewise smooth image prior model and measured CCD camera response
  functions. We also learn the space of noise level functions how noise
  level changes with respect to brightness and use Bayesian MAP
  inference to infer the noise level function from a single image. We
  illustrate the utility of this noise estimation for two algorithms:
  edge detection and feature preserving smoothing through bilateral
  filtering. For a variety of different noise levels, we obtain good
  results for both these algorithms with no user-specified inputs.

Single-Image Noise Level Estimation for Blind Denoising, with Matlab code

  Noise level is an important parameter to many image processing
  applications. For example, the performance of an image denoising
  algorithm can be much degraded due to the poor noise level estimation.
  Most existing denoising algorithms simply assume the noise level is
  known that largely prevents them from practical use. Moreover, even
  with the given true noise level, these denoising algorithms still
  cannot achieve the best performance, especially for scenes with rich
  texture. In this paper, we propose a patch-based noise level
  estimation algorithm and suggest that the noise level parameter should
  be tuned according to the scene complexity. Our approach includes the
  process of selecting low-rank patches without high frequency
  components from a single noisy image. The selection is based on the
  gradients of the patches and their statistics. Then, the noise level
  is estimated from the selected patches using principal component
  analysis. Because the true noise level does not always provide the
  best performance for nonblind denoising algorithms, we further tune
  the noise level parameter for nonblind denoising. Experiments
  demonstrate that both the accuracy and stability are superior to the
  state of the art noise level estimation algorithm for various scenes
  and noise levels.


How many DFT bins must I use to get best precision when analyzing melody scale

You can always improve the accuracy of frequency estimates by quadratically interpolating the peaks in the DFT. This page should explain how to do this: https://ccrma.stanford.edu/~jos/parshl/Peak_Detection_Steps_3.html

You still end up with the issue of converting peaks in the spectrum to pitch, i.e. which peaks correspond to the fundamentals. The fundamental may not even be present as well!


Can anyone suggest good algorithms for CBIR?

A luminosity histogram (especially one that is separated into RGB components) is a reasonable fingerprint for an image - and can be implemented quite efficiently.

Subtracting one histogram from another will produce a new historgram which you can process to decide how similar two images are.
Histograms, because the only evaluate the distribution and occurrence of luminosity/color information handle affine transformations quite well.
Luminosity histograms produce false negatives when the color information in an image is manipulated.


Can someone let me know linearity, time-invariance, causality, memory characteristic of the system?

There's an error in your 3rd section (linearity) - you replace $y(t)$ with $y_1(t)$ on the left side of the equation but not on the right side. You should have
$$y_1(t)=x^2(t-t_0)-\frac{1}{2}\frac{dy_1(t)}{dt}$$.
Now you see that $y_1(t)=y(t-t_0)$.

As for memoryless - start with your initial equation and integrate both sides. Now you have
$$
\int y(t)dt=\int x^2(t)dt-\frac{1}{2} y(t)dt
$$

and rearranging gives
$$
y(t) = 2 \int x^2(t)dt-2 \int y(t)dt
$$
Does that help?


Yes, if you put a pole approaching the unit circle (but not on the unit circle, as this will create an unstable filter) , you will have a bandpass. The center frequency of this bandpass filter depends on how many radians along the unit circle where you have placed the pole (the pole's 'degree').

Moving counter-clockwise from zero will give you an increasing center frequency of the bandpass filter.

Once you start moving towards zero on the Z-plane and away from the unit circle, you will decrease the sharpness Q and consequently increase the bandwidth.


Is there a way in MATLAB to approximate an FIR filter using an IIR Filter?

Sure there is.

The simplest version would be to calculate the transfer function of your FIR filter, sample it on some meaningful frequency grid and then shove it into invfreqz() to design the IIR filter.

Depening on you FIR filter this may work well or be a total disaster. In my experience invfreqz() does often very poorly.

If that's not good enough you can certainly write your own search an optimization algorithm. However, this is NOT a beginner project and requires significant math and coding skills.

See mathworks.com/help/signal/ref/invfreqz.html#mw_3c79f471-99ad-46f8-a680-14c29ff47d07


Filtering signal with added noise (frequencies known)

You want a Wiener filter here.  Basically your frequency response at frequency $f$ should be $\frac{E |S|^2}{E |S|^2 + E |N|^2}$ where $S$ is the Fourier transform of the signal and $N$ of the noise (all dependent on $f$ of course).  So your filters should not aim for a stop band of zero but rather a stop band that will admit an amount of signal based on the expected SNR at that frequency.

There are comparatively few staircases in the total resulting frequency response here, so you might consider going for forward and backward in time IIR filtering (netting a zero phase delay, with normal delay the Wiener filter design becomes pointless for more than a 90degree phase shift) or some FFT with windowed overlap-add filtering.

Wiener filters like that minimize the expected mean square error, so they seem to be the right tool here.


How to test I/Q modulation with only one branch (I or Q)?

The FSW is superheterodyne.

This means that it will not have a DC offset or I/Q imbalance, but the IF path will appear as an asymmetrical lowpass (because it is really a shifted bandpass), which is then compensated for in the calibration.

The first question: have you run the calibration properly?

leave both machines running for half an hour in a stable environment.
disable the generator output (OUTP OFF)
calibrate both (*CAL)

Lots of equipment lenders don't understand that the calibration sequence is used to compensate temperature, run the calibration once (usually on a cold machine) to make the error message go away and tell customers to never touch that button. This gives a fairly bad filter to compensate the IF, and since that filter is applied before downconverting the IF to zero, this will also appear as a complex FIR filter to you.

Second, test your setup with a simple test:

generate a sine wave (*RST, FREQ 2 GHZ, OUTP ON)
look at the vector display in IQ Analyzer (*RST, INST IQ, LAY:REPL '1',VECT)

You should see a single stationary dot there, with a position that corresponds to an arbitrary phase angle that is the sum of the phase between the LOs (they are certainly not the same frequency!) and the digital downconversion inside the FSW, which together give a constant, but only as long as all the oscillators keep running at their fixed frequency.

Once you have reached that state, you should be able to repeat your measurement with more complex signals.


What do troughs represent in a sound wave?

Sound is a pressure wave in a fluid medium. That medium has a "static"or average pressure and sound is simply a small change of that pressure (as a function of time). Negative sound pressure just means that the current pressure is a little lower than the average pressure.
For example, static air pressure is around 100,000 Pascals (Pa) . Normal audible sound is quite a bit smaller, maybe 0.1 Pa (for 76 dB SPL). That means that the pressure are the microphone (or your ear) varies from 99,999.1 Pa to 100,000.1 Pa


Sinogram Fourier Analysis Doesn't Show Clear Bow Tie Shape

Try imtool() and adjust the contrast to show ft similar as in the paper. I set eliminate outliers to 5%, then (contrast) window minimum to 0.0. May I suggest to ask the author about the settings for the display range in his paper.

Now to the signal processing.

The mentioned paper by Vo is public available. Moreover Vo published the algorithm in mathematica code. Note that the code by Vo does a gaussian filtering, i.e. a lowpass of the sinogram. IMHO to achieve similar results as in a real scan.

Generally a real scan shows lowpass behaviour, i.e. the high frequencies in the projections in the sinogram are lost due to the size of the focal point inside the x-ray tube, the detector system and radiation scatter inside the object. In contrast the forward projection algoritm: Those may keep high frequencies in the projections, depending on their base algorithm.

sinogram1 = GaussianFilter[ sinogram1, {{3,2}}];
(* Apply smooth filter: recommend blur more in the vertical direction *)

The Mathematica help tells, that those 3 and 2 steer the xy-sigma of the GaussianFilter(). One can go theroretically through that to find the values of the matrix, but I prefer to run the above filter over a dirac impulse to see the resulting filtermatrix and establish it similar in matlab for use in imfilter(). But I have no access to mathematica. May be s.o can translate the Mathematica GaussianFilter() into a matlab fspecial().

So let's filter the sinogram with a lowpass to ged rid of the high frequencies introduced by radon() :

imtool(ft); %% eliminate 5% outliers and set min to 0.0
h=fspecial('gaussian',[2 3]); %% wider in projection direction
sino_lp = imfilter(sinogram,h,'replicate');
ft_lp=real(fftshift(fft2(sino_lp)));
imtool(ft_lp);

Such a fourier transformed sinogram has lower values in the high frequencies.


What does natural scene image data-set mean?

From one perspective, not much. However, biological perceptual systems have evolved in response to physical properties of natural environments. And such environments possess certain characteristics, such as Laplacian distribution shaped histograms. Therefore natural images are grouped and analyzed together.


Meaning of Phase response of a filter? In simple words?

Here is the same experiment done in my blog: https://poweidsplearningpath.blogspot.com/2020/04/chapter-51-meaning-of-general-linear.html

First, there is a tiny mistake in the question. Not all FIR filters have the properties of Linear Phase. Only the four types of FIR possess the properties. 1 By contrast, all of the IIR filters are not linear phases.

For me, one better description of general linear phase (GLP) is constant group delay. By definition, group delay is negative derivative of phase 2 (In fact, the detail phase/delay relationship can be derived but I suggest we just accept them.). Then, Derivative (group delay) of a linear (phase) is constant and vise versa.

Then, let's move on to the meaning of group delay.
group delay of a frequency represents the delay unit of the filter to that frequency.
So, the filter may treat different frequencies with different delay unit. For an extreme bad example of an non linear filter, the input signal 'do re mi' may becomes 're mi do' in the output. A GLP filter can guarantee that such wield condition will never happen.

Here I wrote an example. The example comes from the Chapter 5.1.2 in the bible of DSP3 and I just implemented the sample.

First, a IIR Filter with phase response like this is given.

Fig. 1.

Here is the group delay (negative derivative of phase response) and the magnitude response. Please note that I delay the frequency in +-0.2 pi for about 150 units. BTW, the filter is a low pass filter so signal higher than 0.8pi is expected to be filtered.

Fig. 2.

Then, let's input a test signal like 'do re mi'. The signal x[n] are '0.8pi, 0.2pi, 0.4pi' in order. The corresponding frequency response is also provided.

Fig. 3.

And here is the output signal. The signal becomes 'empty, 0.4pi, 0.2pi'. The signal component with 0.8pi is filtered out as expected.

Fig. 4.

To make thing more clear, here I points the number of Fig. 3 and Fig. 4 together. For 0.2pi component, the group delay is about 6.39 units but the group delay of 0.2pi component is about 153 unit. The output signal can confirm the prediction from group delay response. That is why the 0.2pi component becomes that last in the output.

Fig. 5.

In summary,

linear phase equals to constant group delay.
GLP FIR filter can guarantee such scenario will never happen. But IIR can never achieve GLP. (But with the same mag frequency spectrm requirement, IIR usually can achieves the spec with lower delay (but not constant) in comparison with FIR.)

Reference:

FIR filter with linear phase, 4 types
https://en.wikipedia.org/wiki/Group_delay_and_phase_delay
A. Oppenheim and R. Schafer, Discrete Time Signal Processing 3rd. 2009

Matlab code

%% System
% H1[z]
b1 = conv([1 -.98*exp(j*.8*pi)],[1 -.98*exp(-j*.8*pi)]);
a1 = conv([1 -.8*exp(j*.4*pi)],[1 -.8*exp(-j*.4*pi)]);
H1 = tf(b1,a1,-1,'Variable','z^-1');

% H2[z]
H2 = tf(1,1,-1,'Variable','z^-1');
for k = 1:4
    ck = 0.95*exp(j*(0.15*pi+0.02*pi*k));
    ck_conj = conj(ck);
    b_tmp = conv([ck_conj -1],[ck -1]);
    b_tmp = conv(b_tmp,b_tmp);
    a_tmp = conv([1 -1*ck],[1 -1*ck_conj]);
    a_tmp = conv(a_tmp,a_tmp);
    H_tmp = tf(b_tmp,a_tmp,-1,'Variable','z^-1');

    H2 = series(H2,H_tmp);
end

% H[z]
H = series(H1,H2);

% Zero-Pole Plot, Fig. 5.2
[b_h,a_h] = tfdata(H );
b_h = cell2mat(b_h);
a_h = cell2mat(a_h);

figure;
zplane(b_h,a_h);
suptitle('Zero-Pole Plot, Fig 5.2');

% System Response.
L=1000;
dw=2*pi/L;
w = -pi:dw:pi-dw;
HH=freqz(b_h,a_h,w);
mag=abs(HH);
phase=angle(HH);

% Fig. 5.3
figure;
subplot(2,1,1);
plot(w,phase);
xticks([-pi -0.8*pi -0.6*pi -0.4*pi -0.2*pi 0 0.2*pi 0.4*pi 0.6*pi 0.8*pi pi]);
xticklabels({'-\pi','-0.8\pi','-0.6\pi','-0.4\pi','-0.2\pi','0','0.2\pi','0.4\pi','0.6\pi','0.8\pi','\pi'});
xlim([-pi pi]);
yticks([-4 -2 -0 2 4]);
ylabel('ARG[H(e^(^j^w^)]');
xlabel('w');
title('Phase response');

subplot(2,1,2);
plot(w,unwrap(phase));
xticks([-pi -0.8*pi -0.6*pi -0.4*pi -0.2*pi 0 0.2*pi 0.4*pi 0.6*pi 0.8*pi pi]);
xticklabels({'-\pi','-0.8\pi','-0.6\pi','-0.4\pi','-0.2\pi','0','0.2\pi','0.4\pi','0.6\pi','0.8\pi','\pi'});
xlim([-pi pi]);
ylabel('arg[H(e^(^j^w^)]');
xlabel('w');
title('Unwrap Phase response');
suptitle('ARG/arg Plot, Fig 5.3');

% Fig. 5.4
figure;
subplot(2,1,1);
plot(w(1:end-1),-1*diff(unwrap(phase))./diff(w));
xticks([-pi -0.8*pi -0.6*pi -0.4*pi -0.2*pi 0 0.2*pi 0.4*pi 0.6*pi 0.8*pi pi]);
xticklabels({'-\pi','-0.8\pi','-0.6\pi','-0.4\pi','-0.2\pi','0','0.2\pi','0.4\pi','0.6\pi','0.8\pi','\pi'});
xlim([-pi pi]);
ylabel('grd[H(e^(^j^w^)]');
title('Group Delay');

subplot(2,1,2);
plot(w,mag);
xticks([-pi -0.8*pi -0.6*pi -0.4*pi -0.2*pi 0 0.2*pi 0.4*pi 0.6*pi 0.8*pi pi]);
xticklabels({'-\pi','-0.8\pi','-0.6\pi','-0.4\pi','-0.2\pi','0','0.2\pi','0.4\pi','0.6\pi','0.8\pi','\pi'});
xlim([-pi pi]);
ylabel('|H(e^(^j^w^)|');
title('Magnitude response');
suptitle('GD/mag Plot, Fig 5.4');

%% Signal
M = 60;
n = 0:M;
w = 0.54-0.46*cos(2*pi*n/M);

N = 512;
x1 = zeros(1,N);
x2 = zeros(1,N);
x3 = zeros(1,N);
dw = 2*pi/N;
w_freq = -pi:dw:pi-dw;

for i = 0:M

    x1(i+M) = w(i+1)*cos(0.2*pi*i);
    x2(i+2*M-1) = w(i+1)*cos(0.4*pi*i-pi/2);
    x3(i+1) = w(i+1)*cos(0.8*pi*i+pi/5);
end
x = x1+x2+x3;
X = abs(fft(x));
X = fftshift(X);

% Fig. 5.5
figure;
subplot(2,1,1);
plot(x);
title('x[n]');
xlim([0,300]);
subplot(2,1,2);
plot(w_freq,X);
xticks([-pi -0.8*pi -0.6*pi -0.4*pi -0.2*pi 0 0.2*pi 0.4*pi 0.6*pi 0.8*pi pi]);
xticklabels({'-\pi','-0.8\pi','-0.6\pi','-0.4\pi','-0.2\pi','0','0.2\pi','0.4\pi','0.6\pi','0.8\pi','\pi'});
xlim([-pi pi]);
ylabel('|H(e^(^j^w^)|');
title('DTFT of X');
suptitle('Input time/Freq., Fig 5.5');

%% Output
y = filter(b_h,a_h,x);

% Fig. 5.6
figure;
plot(y);
xlim([0,300]);
xlabel('n');
title('output y[n], Fig 5.6');

%= Compre the Delay sample point.
figure;
subplot(2,1,1);
plot(x);
xlim([0,300]);
xlabel('n');
ylabel('x[n]');
title('input');
subplot(2,1,2);
plot(y);
xlim([0,300]);
xlabel('n');
ylabel('y[n]');
title('output');


For example a delay of 0.75 sampling periods would have an impulse response like this (red squares sampled from the blue delayed sinc):

Time (horizontal axis) unit is sampling periods. That kind of a filter is not causal as there are non-zero samples at negative times. None of the samples is zero, no matter how far left or right you go. That's what they meant by spreading all over.

About your second question, I have no access to the paper but I know that such all-pass fractional delay filters can be made arbitrarily good (approaching sinc), and the latency increases as you increase filter order. If you want very short latency, you then will also have a very low-order filter, which is of very poor quality as it is completely missing the left side of the sinc. It is a tradeoff between quality and latency, and you probably don't want to trade away all quality.


In a nutshell: yes, it does, at least under certain circumstances.
In a noisy environment, people tend to not only increase their volume when speaking, but also speak in a higher pitch. It's called Lombard Effect, it's well known and you should find plenty of material about it.


Didn't check the code, but the steps you describe in 1. seems correct.
The units of the autocorrelation are time and not frequency, and it is centred around $\tau = 0$. That means that in your $-1024 \leq\tau< 1023$. So you are looking at the wrong peaks.
Now as you can see the autocorrelation doesn't produce a sharp peak which means that it is not a good measure for your case. In general, the autocorrelation doesn't work well with narrow band signals.

I would suggest you try either the MPM or the YIN methods which are based on the autocorrelation, are extremely simple to implement and from my experience yield goo results.


iPhone iOS UIImage how to detect a "laser pointer" dot on a camera feed?

OpenCV is compilable for iOS. While it may not be the most efficient it gives you the option of porting the algorithm. I have done a similar marker tracking problem using the ConDensation algorithm. Look up marker tracking with OpenCV. It is a VERY large research area and the exact algorithm you want depends completely on your application. If I remember correctly there are around 3000 reported image processing techniques - to pick a good one is a real art!

  One solution would be to sample the color of pixels within a certain
  radius of the center and detect a blob bright color. Can anyone
  suggest an algorithm for this activity?

By the way, this is the basic idea behind what is called a particle filter (of which Condensation is one method). Well done you figured out the basic idea by your self!

  Another approach would be to keep track of the average position of the
  dot over the last few frames, thus reducing the guesstimate radius. If
  there's no dot within a pre-defined region, the search region may be
  expanded.

This is called state persistance and can be modelled in various ways. The condensation algorithm uses a stochastic approach not dissimilar to a plain old Kalman filter.

  Is there an efficient way of comparing multiple small images to one
  another to understand if their silhouettes match?

This one is a little more difficult. You could try template matching but I don't know how well it will perform on iOS (very heavy calculation and iOS camera is not well suited to it).


STFT Window Size is 2048, Why is the output 1025?

The continuous Fourier transform possesses symmetries  when computed on real signals (Hermitian symmetry). The discrete version, an FFT (of even length) possesses a slighty twisted symmetry.

The DC coefficient ($F(0)$) is real, as well as the Nyquist one ($F(N/2)$). In between, you get $\frac{2048-2}{2}=1023$ "complex" coefficients, "duplicated" in positive and negative frequencies.

So for real signal, each STFT frame can be represented by $1023+2$ frequency bins, the remaining   1023 being recovered by Hermitian symmetry.

You can get complementary information at FFT of Pure Real Sequences:

  As a result, we can see that if N is even, both $F(0)$ and $F(N/2)$ must
  be real. Given these two values and the complex values $F(1)...F(N/2-1)$,
  (I.E. N numbers in total) the sequence is completely characterised.


FFmpeg audio filter pipeline for speech enhancement

I just came across one article related to the compand command that might be useful (just in case someone looking for help sees this article).

From the article:

  To test this, and my understanding of compand, I added a very simple filter to remove the quiet section of the audio by decreasing their volume dramatically:

ffmpeg -i in.mp3  -filter_complex \
  "compand=attacks=0:points=-30/-900|-20/-20|0/0|20/20" \
  out.wav

  attacks=0 means that I wanted to measure absolute volume, not averaging the sound over a short (or long period of time). When the speaker suddenly yells, or talks back and forth with someone off mic, I want the volume adjustment to be immediate. The downside is that you may hear the volume being clamped.

  points is the actual volume mapping function, and I’ll walk through it:

  I added a mapping of -30/-900, which means that volume below -30db in the original input track gets converted to -900db (completely silent).


What does one need to know in order to build a software synthesizer?

Get Real Sound Synthesis for Interactive Applications, and download the Synthesis Toolkit.  That will give you practical introductions and quality source code for most digital synthesis methods in use.  Then pick whichever method you find most interesting, and start modifying the corresponding STK code to do what you want.  Then try building something similar from scratch.

In addition, like hotpaw2 suggested, anything from JOS is also great introductory material.  In particular, Introduction to Digital Filters with Audio Applications and Mathematics of the Discrete Fourier Transform both start from the basics, and don't assume a lot of background.

From a digital filter/effects perspective, another good starting point is to check out Faust.  You don't even have to learn their block-diagram signal processing language, you can just use the examples to dump out VST, LADSPA, etc. template code as quick way of getting something running [which you can then modify].

That being said, it's even easier to start out with sound rendering, which is basically the same thing, but without worrying about real-time operation or interactivity.  You just save wav files to disk.  It's a great way to experiment with new techniques before going through the trouble of making them run quickly.


Which transform most closely mimics the human auditory system?

In designing such transformations, one should take into account competing interests:

fidelity to the human auditory system (that varies with people), including non-linear or even chaotic aspects (tinnitus)
easiness of the mathematical formulation for the analysis part
possibility to discretize it or allow fast implementations
existence of a  suitable stable inverse

Two recents designs have catch my ears recently:
Auditory-motivated Gammatone wavelet transform, Signal Processing, 2014

  The ability of the continuous wavelet transform (CWT) to provide good
  time and frequency localization has made it a popular tool in
  time–frequency analysis of signals. Wavelets exhibit constant-Q
  property, which is also possessed by the basilar membrane filters in
  the peripheral auditory system. The basilar membrane filters or
  auditory filters are often modeled by a Gammatone function, which
  provides a good approximation to experimentally determined responses.
  The filterbank derived from these filters is referred to as a
  Gammatone filterbank. In general, wavelet analysis can be likened to a
  filterbank analysis and hence the interesting link between standard
  wavelet analysis and Gammatone filterbank. However, the Gammatone
  function does not exactly qualify as a wavelet because its time
  average is not zero. We show how bona fide wavelets can be constructed
  out of Gammatone functions. We analyze properties such as
  admissibility, time-bandwidth product, vanishing moments, which are
  particularly relevant in the context of wavelets. We also show how the
  proposed auditory wavelets are produced as the impulse response of a
  linear, shift-invariant system governed by a linear differential
  equation with constant coefficients. We propose analog circuit
  implementations of the proposed CWT. We also show how the
  Gammatone-derived wavelets can be used for singularity detection and
  time–frequency analysis of transient signals.

The ERBlet transform: An auditory-based time-frequency representation with perfect reconstruction, ICASSP 2013

  This paper describes a method for obtaining a perceptually motivated
  and perfectly invertible time-frequency representation of a sound
  signal. Based on frame theory and the recent non-stationary Gabor
  transform, a linear representation with resolution evolving across
  frequency is formulated and implemented as a non-uniform filterbank.
  To match the human auditory time-frequency resolution, the transform
  uses Gaussian windows equidistantly spaced on the psychoacoustic “ERB”
  frequency scale. Additionally, the transform features adaptable
  resolution and redundancy. Simulations showed that perfect
  reconstruction can be achieved using fast iterative methods and
  preconditioning even using one filter per ERB and a very low
  redundancy (1.08). Comparison with a linear gammatone filterbank
  showed that the ERBlet approximates well the auditory time-frequency
  resolution.

And I shall mention also:

An Auditory-Based Transform For Audio Signal Processing, WASPAA 2009

  An auditory-based transform is presented in this paper. Through an
  analysis process, the transform coverts time-domain signals into a set
  of filter bank output. The frequency responses and distributions of
  the filter bank are similar to those in the basilar membrane of the
  cochlea. Signal processing can be conducted in the decomposed signal
  domain. Through a synthesis process, the decomposed signals can be
  synthesized back to the original signal through a simple computation.
  Also, fast algorithms for discrete-time signals are presented for both
  the forward and inverse transforms. The transform has been approved in
  theory and validated in experiments. An example on noise reduction
  application is presented. The proposed transform is robust to
  background and computational noises and is free from pitch harmonics.
  The derived fast algorithm can also be used to compute continuous
  wavelet transform


Analog Filter Prototype to Direct Form Second Order Cascade

The coefficients in the function lowpass3() are the numerator and denominator coefficients of the two sections (first and second order) of the normalized continuous-time transfer function. The continuous-time transfer function looks like this:

$$H(s)=\frac{b_{21}s^2+b_{11}s+b_{01}}{s^2+a_{11}s+a_{01}}\cdot\frac{1}{s+a_{02}}\tag{1}$$

Note that this normalized filter has  a cut-off frequency of $1$. The two functions tf2s and tf1s transform these two sections to the discrete-time domain using the bilinear transform, whereby the normalized cut-off frequency of the continuous-time filter is mapped to the desired cut-off frequency in the discrete-time domain. So whereas the coefficients in $(1)$ are independent of the cut-off frequency (because the filter is normalized), the coefficients of both final filter sections in the discrete-time domain do depend on the cut-off frequency.

Your example of using ncauer followed by zp2sos will not give you the coefficients of the discrete-time filter with the desired cut-off frequency. Instead, you just the compute the coefficients of the first and second order sections of the normalized continuous-time low pass filter, i.e., the coefficients hard-coded in lowpass3(). What you miss in that case is the transformation step from continuous-time to discrete-time, and the mapping of the normalized cut-off frequency to the desired cut-off frequency.


http://www.mathworks.com/help/toolbox/signal/ref/spectrogram.html

This pretty standard fare. For audio I would recommend a hanning() window and an overlap of half a window length so

x = randn(44100,1); % some random noise sample of a second
S = spectrogram(x,hanning(512),256); % short term FFT

Assuming a sample rate of 44.1 kHz this gives you frames about 11.6 ms length. You can get to exactly 10ms but that requires FFTs that are not a power of 2.


Harmonic distortion of sine signal. plotting resulting spectrum

Here is a way that you can approach this mathematically.  It's not the only way.  Another way to approach this is with Tchebyshev polynomials and I will leave that to someone else.

This approach requires understanding two ingredients.

an identity that is a consequence of Euler's formula:

$$ \cos(\omega t) = \tfrac12(e^{j\omega t}+e^{-j\omega t}) $$
or
$$ \sin(\omega t) = \tfrac1{2j}(e^{j\omega t}-e^{-j\omega t}) $$

the binomial theorem:

$$\begin{align}
 (a+b)^N &= \sum_{n=0}^N {N \choose n} \, a^{N-n} \, b^n \\
         &= \sum_{n=0}^N {N \choose n} \, a^n \, b^{N-n} \\
\end{align}$$
or
$$\begin{align}
 (a+b)^N &= \sum_{n=0}^N \frac{N!}{(N-n)! \, n!} \, a^{N-n} \, b^n \\
         &= \sum_{n=0}^N \frac{N!}{(N-n)! \, n!} \, a^n \, b^{N-n} \\
\end{align}$$

Now you might see where this is going.  So using the cosine identity (so i don't have to mess around with powers of $j$ in the sine identity):

$$\begin{align}
 \cos^k(\omega t) &= \frac{1}{2^k}(e^{j\omega t}+e^{-j\omega t})^k \\
\\
 &= \frac{1}{2^k} \sum_{n=0}^k \frac{k!}{(k-n)! \, n!} \, (e^{j\omega t})^{k-n} \, (e^{-j\omega t})^n \\
\\
 &= \frac{1}{2^k} \sum_{n=0}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-n)\omega t} \, e^{-jn\omega t} \\
\\
 &= \frac{1}{2^k} \sum_{n=0}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} \\
\end{align}$$

It appears that we will have to deal with the two cases of even $k$ vs. odd $k$ differently.

For even $k \ge 2$:

$$\begin{align}
 \cos^k(\omega t) &= \frac{1}{2^k}(e^{j\omega t}+e^{-j\omega t})^k \\
\\
 &= \frac{1}{2^k} \sum_{n=0}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} \\
\\
 &= \frac{1}{2^k} \left( \sum_{n=0}^{\tfrac{k}{2}-1} \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} + \frac{k!}{\tfrac{k}{2}! \, \tfrac{k}{2}!} + \sum_{n=\tfrac{k}{2}+1}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} \right) \\
\\
 &= \frac{1}{2^k} \left(\frac{k!}{\left(\tfrac{k}{2}!\right)^2} +  \sum_{n=0}^{\tfrac{k}{2}-1} \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} + \sum_{n=0}^{\tfrac{k}{2}-1} \frac{k!}{n! \, (k-n)!} \, e^{j(k-2(k-n))\omega t} \right) \\
\\
 &= \frac{1}{2^k} \frac{k!}{\left(\tfrac{k}{2}!\right)^2} + \frac{1}{2^k} \sum_{n=0}^{\tfrac{k}{2}-1} \frac{k!}{(k-n)! \, n!} \, \Big( e^{j(k-2n)\omega t} + e^{-j(k-2n)\omega t} \Big) \\
\\
 &= \frac{1}{2^k} \frac{k!}{\left(\tfrac{k}{2}!\right)^2} + \frac{1}{2^k} \sum_{n=0}^{\tfrac{k}{2}-1} \frac{k!}{(k-n)! \, n!} \, 2 \cos\big((k-2n)\omega t\big) \\
\\
 &= \frac{1}{2^k} \frac{k!}{\left(\tfrac{k}{2}!\right)^2} + \frac{2 k!}{2^k} \sum_{n=1}^{\tfrac{k}{2}} \frac{1}{\left(\tfrac{k}{2}+n\right)! \, \left(\tfrac{k}{2}-n\right)!} \, \cos(2n\omega t) \\
\end{align}$$

Note that only DC and even-numbered harmonics come outa the even powers of $\cos^k(\omega t)$.

For odd $k \ge 1$:

$$\begin{align}
 \cos^k(\omega t) &= \frac{1}{2^k}(e^{j\omega t}+e^{-j\omega t})^k \\
\\
 &= \frac{1}{2^k} \sum_{n=0}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} \\
\\
 &= \frac{1}{2^k} \left( \sum_{n=0}^{\tfrac{k-1}{2}} \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} + \sum_{n=\tfrac{k+1}{2}}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} \right) \\
\\
 &= \frac{1}{2^k} \left( \sum_{n=0}^{\tfrac{k-1}{2}} \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} + \sum_{n=0}^{\tfrac{k-1}{2}} \frac{k!}{n! \, (k-n)!} \, e^{j(k-2(k-n))\omega t} \right) \\
\\
 &= \frac{1}{2^k} \sum_{n=0}^{\tfrac{k-1}{2}} \frac{k!}{(k-n)! \, n!} \, \Big( e^{j(k-2n)\omega t} + e^{-j(k-2n)\omega t} \Big) \\
\\
 &= \frac{1}{2^k} \sum_{n=0}^{\tfrac{k-1}{2}} \frac{k!}{(k-n)! \, n!} \, 2 \cos\big((k-2n)\omega t\big) \\
\\
 &= \frac{2 k!}{2^k} \sum_{n=0}^{\tfrac{k-1}{2}} \frac{1}{\left(\tfrac{k+1}{2}+n\right)! \, \left(\tfrac{k-1}{2}-n\right)!} \, \cos\big((2n+1)\omega t\big) \\
\end{align}$$

Note that only odd-numbered harmonics come outa the odd powers of $\cos^k(\omega t)$.

So now what if we pass this sinusoid:

$$ x(t) = A \, \cos(\omega t) $$

through this polynomial mapping function:

$$\begin{align}
 y(t) &= a_0 + a_1 x(t) + a_2 x^2(t) + a_3 x^3(t) + ... + a_{K-1} x^{K-1}(t) + a_K x^K(t) \\
      &= \sum_{k=0}^{K} a_k \big( x(t) \big)^k \\
\end{align} $$

?

I believe we have to split the even-order terms from the odd-order terms in the summation

$$\begin{align}
 y(t) &= \sum_{k=0}^{K} a_k \big( x(t) \big)^k \\
\\
  &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \big( x(t) \big)^{2k} + \sum_{k=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} \big( x(t) \big)^{2k+1} \\
\\
  &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \big( A \cos(\omega t) \big)^{2k} + \sum_{k=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} \big( A \cos(\omega t) \big)^{2k+1} \\
\\
  &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} A^{2k} \cos^{2k}(\omega t) + \sum_{k=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} A^{2k+1} \cos^{2k+1}(\omega t) \\
\\
  &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} A^{2k} \left( \frac{1}{2^{2k}} \frac{(2k)!}{(k!)^2} + \frac{2 (2k)!}{2^{2k}} \sum_{n=1}^{k} \frac{1}{\left(k+n\right)! \, \left(k-n\right)!} \, \cos(2n\omega t) \right) + \sum_{k=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} A^{2k+1} \left( \frac{2 (2k+1)!}{2^{2k+1}} \sum_{n=0}^{k} \frac{1}{(k+1+n)! \, (k-n)!} \, \cos\big((2n+1)\omega t\big) \right) \\
\\
  &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{A^{2k}}{2^{2k}} \frac{(2k)!}{(k!)^2} + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} \sum_{n=1}^{k} a_{2k} \frac{2 A^{2k} (2k)!}{2^{2k}\left(k+n\right)! \left(k-n\right)!} \, \cos(2n\omega t) \qquad + \sum_{k=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} \sum_{n=0}^{k} a_{2k+1} \frac{2 A^{2k+1} (2k+1)!}{2^{2k+1}(k+1+n)!(k-n)!} \, \cos\big((2n+1)\omega t\big) \\
\\
  &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{A^{2k}}{2^{2k}} \frac{(2k)!}{(k!)^2} + \sum_{n=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} \sum_{k=n}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{2 A^{2k} (2k)!}{2^{2k}\left(k+n\right)! \left(k-n\right)!} \, \cos(2n\omega t) \qquad + \sum_{n=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} \sum_{k=n}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} \frac{2 A^{2k+1} (2k+1)!}{2^{2k+1}(k+1+n)!(k-n)!} \, \cos\big((2n+1)\omega t\big) \\
\end{align} $$

$\lfloor u \rfloor$ is the floor() function, which is the largest integer no larger than the argument $u$.

$$ \lfloor u \rfloor \le u  < \lfloor u \rfloor + 1 $$
or
$$  u-1 < \lfloor u \rfloor \le u $$

Also note that if $K$ is even
$$ \left\lfloor \tfrac{K-1}{2} \right\rfloor < \left\lfloor \tfrac{K}{2} \right\rfloor = \tfrac{K}{2} \quad \text{and} \quad \left\lfloor \tfrac{K-1}{2} \right\rfloor < \tfrac{K-1}{2} $$

If $K$ is odd
$$ \left\lfloor \tfrac{K}{2} \right\rfloor = \left\lfloor \tfrac{K-1}{2} \right\rfloor = \tfrac{K-1}{2} \quad \text{and} \quad \left\lfloor \tfrac{K}{2} \right\rfloor < \tfrac{K}{2} $$

$$\begin{align}
 y(t) &= \sum_{k=0}^{K} a_k \big( A \cos(\omega t) \big)^k \\
\\
     &= \sum_{n=0}^{K} b_n \cos(n \omega t) \\
\end{align}$$

where, apparently the DC component of the Fourier cosine series is:

$$ b_0 =  a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{(2k)!}{(k!)^2} \left(\frac{A}{2}\right)^{2k}  $$

The coefficient for the even harmonic terms:

$$ b_{2n} = 2 \sum_{k=n}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{(2k)!}{\left(k+n\right)! \left(k-n\right)!} \left(\frac{A}{2}\right)^{2k}  $$

And the coefficient for the odd harmonic terms:

$$ b_{2n+1} = 2 \sum_{k=n}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} \frac{(2k+1)!}{(k+1+n)!(k-n)!} \left(\frac{A}{2}\right)^{2k+1}  $$


What is the most efficient type of software bandpass filter?

Frequency domain filtering (FFT), as suggested by some comments, is definitely wrong -- it's even slower, or same speed at best! A recursive filter (IIR) is the fastes possible solution. If you choose a typical second order filter (called biquad in engineering slang) of Butterworth type and do your math right (factoring out coefficients) you only have 3 multiplications and five addition.

Edit: link to classics page: filter cookbook


Get difference equation of a filter from MATLAB transfer function coefficients

There are two problems with your code. First, you should not invert every other $b$ coefficient, they should all have a positive sign, just like the $a$ coefficients that all have a negative sign. Second, - and that's the main problem -, why do you quantize the coefficients? Your actual coefficient $b[1]$ is 0.00266475348894505 (and even longer), and you rounded it to 0.0027. That's bad for the $b$ (numerator) coefficients, but detrimental for the $a$ (denominator) coefficients, because the filter's characteristic is significantly changed, and the filter can become unstable. You should just use the coefficients as stored in the vectors A and B.

Implementing the difference equation is a nice exercise to understand what's going on, but I think you'd be better off using filter.m.


Why is being jointly WSS important in signal estimation with LMMSE estimator?

Recall that in the derivation of the linear MMSE estimator we force the error term $\hat{s}_k-s_k$ to be orthogonal to the signal $y_k$ i.e. $E[(\hat{s}_k-s_k)y_{k-j}] = 0$ for all $j$. The assumption of joint wide sense stationarity allows one to write $E[s_k y_{k-j}]$ as a function of the lag $j$, i.e. the covariance isn't a function of time $k$, just the lag. As a result, the optimal filter coefficients $h_k$ can be obtained by solving a system of linear equations.


Yes, it's perfectly possible.  All that you'll need is to model how you think the angular velocity components of the state will evolve.  Usually simple brownian (random) motion is enough, at least to start with.  If you know more about how the angular velocities are constrained, then you can include that in the model.

All that it means is there are no outputs from your signal model directly from those elements of the state.

For example, in this answer to another question the state is made up of the $x$ position and the $x$ velocity, but there is no velocity measurement.

$$
\begin{align*}
\mathbf{x}_{k}  & =\left(\begin{array}[c]{c}x_{k}\\
\dot{x}_{k}\end{array}
\right)
\end{align*}
$$

All that means is the output matrix $H$ zeros out that part of the state when producing the measurement.


Properties of Discrete-time Sinusoidal Signal

Any discrete-time periodic signal with period $N$ can be written in terms of its DFT coefficients $X_k$:

$$x(n)=\frac{1}{N}\sum_{k=0}^{N-1}X_ke^{j(2\pi/N) nk}$$

So whatever is true for a single sinusoid (or complex exponential) remains true for a general periodic signal.

Two discrete-time sinusoids with a frequency separation of $2\pi k$ are identical because the sine function is periodic with period $2\pi$:

$$x(n)=\sin(\omega n)\\
y(n)=\sin[(\omega+2\pi k)n]=\sin(\omega n +2\pi k n)=\sin(\omega n)=x(n)$$

Finally, the maximum frequency for a discrete-time sinusoid (and consequently for any discrete-time signal) is $\omega=\pi$:

$$\cos(\pi n)=(-1)^n$$

If the frequency is further increased you will observe aliasing:

$$\cos[(\pi+\Delta\omega)n]=\cos[(\pi+\Delta\omega)n-2\pi n]=\cos[(\Delta\omega-\pi)n]=\\=\cos[(\pi-\Delta\omega)n],\quad 0\le\Delta\omega\le\pi$$

So from the maximum frequency of $\pi$ you get a decrease in frequency down to frequency $0$ (for $\Delta\omega=\pi$).


What input to use for an AR model of a vowel sound?

A pitch estimator is commonly used to find the vocal periodicity.  Common pitch estimators include cepstrum/cepstral analysis, harmonic product spectrum, and composite algorithms, such as YAAPT.


Autocorrelation function of a triangular wave

I couldn't quite follow your code (you calculate two different versions of the ACF?), but I believe the problem with the plot being shifted toward zero is that xcov calculates the cross covariance not correlation. Remember that cross covariance is subtracting the mean, so there should be a shift! You should be calling xcorr instead.
The following code snippet reproduces the Figure 9 from your screenshot.
tau = 1;
b = 100;
Em = sqrt(6);
T1 = 2*b;
signal = zeros(T1, 1);
step = Em/(b+1);
signal(tau:tau+b) = linspace(0, Em, b+1);
signal = repmat(signal, 1000, 1);
[xcorrOut, lag] = xcorr(signal, 'unbiased');

figure
hold on
grid on
plot(lag, xcorrOut)
xlim([-5*b, 5*b])

Edit
In response to the new edit in the question, the OP wants a way to convert the cross correlation result to a cross covariance result. The cross covariance can be obtained from the cross correlation by subtracting the product of the mean and conjugate mean, https://en.wikipedia.org/wiki/Cross-covariance:
$$ \varphi_{11}(\tau) - |\mu_{f_1}|^2$$
where $\mu_{f_1}$ is the mean of $f_1(t)$, as shown in Figure 7 of the screen capture. This additional code snippet builds upon the above to show how adjusting the correlation can get you the same result as computing the covariance.
mu = mean(signal);
xcovFromXcorr = xcorrOut - mu*conj(mu);
[xcovOut, lag] = xcov(signal, 'unbiased');

figure
hold on
grid on
plot(lag, xcovFromXcorr, 'LineWidth', 2)
plot(lag, xcovOut, '--', 'LineWidth', 2)
xlim([-5*b, 5*b])


Timing error detector used for symbol timing synchronization

It depends on your system :)

Seriously, it depends. Gardner is a pretty solid algorithm with a simple implementation, especially if you're using BPSK or QPSK. It can even be used in some MIMO systems. Its main disadvantage is that it requires two samples per symbol. M&M requires only one, but it has its own disadvantages.

I think that these days the preferred solution is some variation of the polyphase detectors first proposed by fred harris. However, if you're experimenting and building your own system and you want to keep it simple, Gardner is a safe bet.


How to interpret the window size and FFT size while performing FFT?

There are lots and lots  of related questions on this topic in this SE and also on many discrete signal processing books.

https://www.sjsu.edu/people/burford.furman/docs/me120/FFT_tutorial_NI.pdf
https://holometer.fnal.gov/GH_FFT.pdf


HINT:

I assume that $n$ is an integer and that we're talking about discrete-time signals. You need to simplify the expressions for $h(n)$ by noting that $\delta(n)=1$ for $n=0$ and zero otherwise. So, for which values of $n$ does $\sin(1+|n|)$ equal zero, and for which values of $n$ does $\sin(2\pi n/N)$ equal zero?


How to determine the step response given a transfer function?

The short answer is that the step response is the inverse Laplace Transform of $\frac{1}{s}G(s)$.

Here are some hints that should help in being able to solve the problem. Since it is likely a homework problem I wanted to provide more details as relevant background and steps toward a solution rather than detail the mathematics for the specific solution. A plot of the resulting step response is included at the end to validate the solution.

First, a review of the simple case of real negative poles on the s-plane along with their inverse (unilateral) Laplace Transform.

As shown, a pole given by the transfer function $H(s) = \frac{1}{s+\alpha}$ has an inverse Laplace transform $x(t) = e^{-\alpha t}$ for $t\ge 0$. This is easy to prove but is also well tabulated in all tables of the Laplace Transforms such as at the bottom of this page: https://en.wikipedia.org/wiki/Laplace_transform

Specifically if $H(s)$ describes the transfer function of the system (output divided by input when the units are $s$), $h(t)$ is the impulse response of the system. As an aside demonstrating one utility of the Laplace Transform, if we restrict $s$ to just be the imaginary axis, we get the frequency response $H(j\omega)$ which in this case would be a low pass filter. Whether $s$ is real or imaginary, we typically refer to the s-domain as the "frequency domain" and all points on the surface are given as pole locations by the transforms of $e^{st}$. So with $h(t)$ being the impulse response, what we see as graphed is the response in time at the output if an impulse was presented at the input. Notice that the pole closest to the imaginary axis takes the longest to decay in its impulse response. This is a key take-away and is the reason we call the negative poles closest to the axis "dominant poles". Like an unwanted guest, they take way too long to go away. Poles further away from the imaginary axis decay much more rapidly in time.

Now consider the summation of these two transfer functions:

$$G(s) = H_1(s) + H_2(s)$$

A summation of two systems in the frequency domain is a parallel process in their signal flow implementation as depicted in the figure below. In contrast the function $H_1(s)H_2(s)$ would be implemented as the cascade of the two systems in series.

Now for the system above represented by the summation of our two systems, the total output due to an impulse at the input would be the sum of the impulse response from each system. Here we see how the pole that was further away from the imaginary axis quickly decays so does not significantly contribute to the final settling time of the output. This settling time (and the overall shape of the response if the poles are far enough apart) is dominated by the pole closest to the imaginary axis, our dominant pole.

Further this combined system is given as:

$$G(s) = \frac{1}{s+b} + \frac{1}{s+a} = \frac{2s+a+b}{(s+a)(s+b)}$$

Notice how adding them together introduced a zero (at $z_1=-(a+b)/2$) but did not change the pole locations. Importantly this shows how given the general combined expression on the right side of the above equation, the process of Partial Fraction Expansions can be used to put this in the form of the parallel (summed) system as described by the left side of the above equation. This is then very easy to solve for the inverse transform,  providing the impulse response, here given as:

$$g(t) = e^{-bt} + e^{-at}$$

Step Response

The above shows how to solve for the impulse response, which is the inverse Laplace Transform of the transfer function $G(s)$.  To get the step response note that the Laplace Transform of an integrator is $1/s$ (also in the Wikipedia link above and a very important takeaway). Similarly the Laplace Transform of a differentiator is $s$ (less the initial condition). Yes commit these to memory as they will come up again and again and again. So imagine if we prefaced our system under test with an integrator: The integratiion of an impulse is a step! So for any given system, if we simply multiply it's transfer function by $1/s$ (which means putting an integrator in cascade or series with the system), the output defined by the inverse Laplace Transform of that result will be the step response! It's that simple. Taking that further if we multiplied by $1/s^2$ we would get a ramp response, etc.

Solving this again is aided by Partial Fraction Expansion:

$$\frac{1}{s}H_1(s) = \frac{1}{s(s+b)}= \frac{1}{b}\bigg(\frac{1}{s} -\frac{1}{s+b} \bigg)$$

From which we can easily derive the inverse Laplace Transform providing the step response $s(t)$:

$$s(t) = \frac{1}{b}\mathscr{L}^{-1}\bigg(\frac{1}{s} -\frac{1}{s+b} \bigg) = \frac{1}{b}(1-e^{-bt})$$

To easily see how the consideration of the dominant pole still applies, consider that a step at the input to the system is an infinite series of impulses each at time $t$ (and to be accurate each is weighted by $dt$ but that is just a scaling). Each impulse produces it's related impulse response, with the responses from llater impulses (accumulating) with the tails of the responses from previous impulses. The accumulation of the weighted impulse responses is described mathematically as convolution but results in the step response at the output. Therefore the dominant pole would also dominate the resulting step response.

With that background in mind, here are the main points specific to solving the OP's problem:

The step response can reasonably estimated from the dominant pole alone.

For an accurate and complete solution, use Partial Fraction Expansion to factor decompose $G(s)$ into three parallel systems $G_1(s)$, $G_2(s)$, $G_3(s)$, which can then each be easily solved for their inverse Laplace Transforms (this would provide the impulse response of $G(s)$.

The step response for each can be solved as described above by first multiplying by $1/s$ and then solving for the inverse Laplace Transform. The result for this is plotted below with $G_1(s)$ containing the dominant pole $s=-0.1$, $G_2(s)$ containing the weaker pole  at $s=-100$ and $G_3(s)$ containing the two repeated poles at $s=-100$ (consistent with Partial Fraction Expansion for repeated poles).

It can be confirmed how closely $G_1(s)$ follows the step response for the system with only the dominant pole given by $1/(s+0.1)$. Also notice the significant difference in time and amplitude scales between the three plots. $G_1(s)$ takes over 50 seconds to settle while $G_2(s)$ and $G_3(s)$ are nearly settled at 50 milliseconds. The sum of all three of these will closely match $G_1(s)$ alone. Further the final values for $G_2(s)$ and $G_3(s)$ are significantly smaller than the final value for $G_1(s)$, but the sum of all three final values will match the final value for combined system $G(s)$.

Comparing the overall step response to that given by the dominant pole highlights the significance of recognizing the dominant pole and how that can be used to simplify most practical problems (I have time now to write this all out from all the time I saved before by recognizing the dominant pole!):

The more significant difference between the two is in the very beginning, so zooming in to show that detail:

The MATLAB/Octave with the control toolbox to plot the step response for $G(s)$ is:

num = -[1 -3]*10000/3;
den = conv([1 0.1], [1 200 10000]);
G = tf(num, den)
step(G);


Is the basis of the sparse signal assumed known in compressed sensing?

Compressive Sensing is an approach to reconstruct sparse signals from incomplete set of measurements.

  In doing so, we need to know $Ψ$ to recover $x$, don't we?

yes, we do.

  But if we know $Ψ$, why don't we just measure
  $z=Ψ^{−1}x=Ψ^{−1}Ψα=a$, and store the value and index of the
  non-zero components of $z (=α)$?

Compressive Sensing  as the name indicates is a method to sense the incoming signal not to compress it. Which means the signal is not present yet and we are going to sample it, this sampling is efficient and requires minimum number of samples for signal reconstruction. There is a misunderstanding here, we might know $Ψ$ in advance but we usually do not have the signal itself $x$, if we do, clearly compression techniques are preferred. The other reason is that sampling domain could be different from the domain in which the signal has a sparse representation. For example, we might sample a fourier sparse signal in time (e.g. radiowaves) or sample a fourier sparse image in space (e.g. Cameras, MRI,..). And then reconstruct it it in fourier domain.

  Can the dimensions ($m$) of $y$ be less than $2k$?

Observing this bands guarantees perfect reconstruction (with overwhelming high probability) for any $k$-sparse signal, but reconstruction is still probable for $m < 2k $; So, as the number of measurements ($m$) is reduced probability of reconstruction reduces as well. But if you mean maintaining the perfect reconstruction probability while minimizing $m$, I think it would depend on your matrix design.


Why use a gaussian low pass filter, when we can convolve with gaussian kernel

In practice gaussian low-pass filters are written using F.I.R. or I.I.R. filters, first the filter runs left then right on each horizontal line, then the filter runs up then down on each vertical line.  It's much cheaper than doing FFTs and IFFTs.

Your two operations are not performing the same operation.  The effects of the filtering around the edges will be very different.  You are doing a circular convolution.  For the Discrete Fourier Transform the convolution theorem says that the multiplication of fourier transforms is equivalent to convolving the image infinitely repeated/tiled in the plane.  So around the left edge, for example you are blurring together some components from the right side of the picture.  (For example try performing your operations on a picture that is black on the left half and white on the right half.)

Your "low pass" filter is padding the image and the kernel to try to reduce some of the edge effects.  It's more expensive, but it is closer to what you would get if you did a "real" convolution of a gaussian with your image.

In your application the circular convolution effect of blurring in data from the other side of the image may be okay, but in many application it is unacceptable.

As to the reason that the FFT of an $N\times N$ image is 10 times faster than the FFT of a $2N \times 2N$ image.  A 2d FFT is done by taking 1d FFTs of each horizontal line, then 1d FFTs of each vertical line.  So in the $N \times N$ case you are doing $2N$ 1-d FFTs, each of which take $O(N \log N)$ for a total of $O(2N^2 \log N)$.  In the $2 N \times 2 N$ case you are doing $4N$ 1-d FFTs each of which take $O(2N \log 2N)$ for a total of (a little more than) $O(8N^2 \log N)$.  Add in cache effects for the larger image and 10x slower seems quite reasonable.


What is the effect of quantization bucket width in Locality Sensitive Hashing?

I think you've got the point, but that you are getting confused by the additional tables that you think should be used to implement buckets.

Let's rephrase the quote from the tutorial first:

  In LSH, you create a partition of the entire space. If you use a large
  bucket width $w$, then obviously you need less buckets to span the
  entire space and thus you get shorter codes. Since you have less
  buckets, your hash table has less slots (obvious) but there are also
  more points per bucket.

Note that by definition of LSH all the points $x$ that fall into a given bucket $B$ have the same binary code:
$$\forall x \in B, \text{LSH}(x) = b,$$
where $b$ is some constant.

Now, let's consider what happens when you make a query on the hash table: if you use $b$ as query to the hash table, you will get all the points $x$ in the above formula.
The only information that you have on these points is that their binary code is equal to $b$, you do not have additional information: they may be all around the bucket, or in its center, etc. So, you can't avoid comparing your exact query point with each of them (using their real coordinates) if you want to find a special neighbour (usually the closest one).

Since you need to loop over all the points inside the bucket, this search is linear, and here comes the trade-off mentioned in the tutorial.

When you implement some LSH, you are free to use whatever technique you want for both the hash table and the buckets.
In my hashing implementation, the hash table was actually a table (for simplicity, because finding a bucket becomes then a simple table lookup), but the buckets themselves were implemented as lists for memory efficiency (you allocate a slot only when you add a point to the bucket).

Your question implies that you would implement the buckets themselves as arrays, but this may not always be desirable when you do not know anything about the distribution of the points in your data set or when your programming framework (in my case, vanilla C/C++) does not allow you to easily resize arrays.


Limitation on the shift theorem of DFT due to frequency resolution?

But what if I want to shift different frequency components by a different number of samples? In that case, you can't find FFT of the whole data set and multiply by the same exponential. You would have to treat different frequencies separately right?

What you're describing is a system with a non-constant (i.e. frequency-dependent) group delay.

Group delay is the derivative of phase, hence,  $\frac{\mathrm d\phantom f}{\mathrm df}\phi(f)$ of a system.

In filter design, you (most of the time) want a constant group delay, so you design your filter with linear phase (which means you need to use FIR filters, and you need to make them symmetric); the derivative of a linear function is a constant!

Now, you want a given "group delay profile". Sooo, what about actually implementing the following idea:

Essentially, you approach this via the window design method that one typically uses to design linear-phase filters of a wanted amplitude response. You, instead, want a to define a group delay, and hence, a phase response over frequency.

Convert a definition of your group delay to a phase definition by numerically integrating that profile (i.e. do a cumulative sum). (the integration is the inverse operation to the differentiation that delivers the group delay from the phase response; so, the result should be that each value is a $e^{j2\varphi}$, with $2\pi\varphi$ being the relative delay you want). Do that definition at a high "oversampling", i.e. if you just need an actual "resolution" for your frequency-dependent delay of sampling rate / 10, define that function at 200 points. Note down that you have a "frequency-sampling-rate" of $r_\text{freq}=\frac{f_\text{sample}}{200}$.
Convolve the result with a low-pass linear filter impulse response. Why will become clearer in the next step.
Convert the result to time domain by doing an IDFT (IFFT). The effect of the low-pass filter above is now, a lot of the resulting values are (close enough to) zero. So, you can omit these.

If I'm not mistaken, you now have a filter of desired group delay profile plus a constant (the constant stems from the low-pass filtering in frequency domain).

Another route if you only care about non-overlapping bands and these should have very different group delays, making it hard to properly build that filter via the filter design method:

Find the "hardest" (i.e. narrowest) bands; design a FIR band pass filter for these using your favorite math toolkit's filter design methodology (I propose designing an equivalent-width low pass and shifting it in frequency). The length of that filter will be the longest you'll need for all the band pass filters. Sufficiently over-sample, i.e. re-design that filter, but use let's say 20 times the length. Use that length for all other filters, too.
The group delay of a linear-phase phase filter is half its length. Now, all your bandpasses have the same group delay $\tau_\text{offset}$.
Adjust the group delay of the individual filters by transforming (FFT), multiplying by a complex sinusoid of appropriate frequency, and transforming back (IFFT). If your frequency isn't (close enough) to $f_s / \text{length}(filter)$, you'll get leakage / Gibb's phenomenon. You can trade accuracy for oscillations by applying a window. If your filters are long, it's probably not even necessary to window.
Sum up all the filter impulse responses.
(optional) if it turns out that a lot of the resulting taps are close to zero, truncate.


Is is necessary to learn digital signal processing before starting research in digital image restoration?

My advice would be: do not delay the start of your research, but integrate the study of signal processing basics in the time you have for your research. Basic signal processing knowledge is essential for a proper understanding of the methods of digital image restoration. Divide your time between learning any basic knowledge you might be missing, and studying specific literature on image restoration.


applicability of dictionary learning in compressed sensing problem

Why of course, it is certainly possible. As you already noted, dictionary learning is tightly coupled to sparse recovery, which is what you need for compressed sensing.

In DL, we are trying to factor $Y = D X$, which is only possible if additional assumptions are made. Typically that includes some notion of sparsity of $X$ and some additional normalization constraints on $D$. There are a few more conditions though to even render the problem unique. Note that $D$ might not be a basis, in general it is not even square. Typical strategies for DL include alternating optimization (have a look for MOD-based algorithms): starting from some guess about one of them, estimate the other one, and then take turns. This takes advantage of the fact that while finding both $D$ and $X$ is a hard (non-convex) problem, finding one of them given the other is much easier:

Finding $D$ given $X$ (and $Y$) is a linear least squares problem
Finding $X$ given $D$ (and $Y$) is a sparse recovery problem (which by itself is also non-convex but has nice convex relaxations that are quite well studied).

In the CS setup, we are encountering a similar situation, except that here, $D$ is typically composed of a measurement matrix we know (say, $\Phi$) and some sparsifying dictionary or basis (say, $A$) we may or may not know so that $D=\Phi A$ and $Y=D X = \Phi A X$. In this case $A$ could be a basis (i.e., it is square) but it could also be some more general dictionary (i.e., flat). The matrix $\Phi$ on the other hand is always flat. Given that, the above strategy still applies and you can apply alternating optimization:

Finding $A$ given $\Phi$, $X$, and $Y$ is a linear least squares problem (unless you need $A$ to be a basis, which adds an orthonormality constraint)
Finding $X$ given $\Phi$, $A$, and $Y$ is a sparse recovery problem (in no way different from what is already used in the CS context).

That said, there are strategies to optimize $D$ and $X$ jointly, like the K-SVD. They might not be directly applicable.


machine learning on multiple signal at once

You say you're "a beginner" so I will suggest using Keras (https://keras.io/). It is written in Python and runs on top of Tensorflow which is a neural network toolbox written by Google. Keras allows you to build a network layer-by-layer in an easy way. For example, here is an example of building a 1D CNN (https://keras.io/examples/imdb_cnn/), please note that the example contains many things that you do not need or are optional. My advice, start with a simple network and know how to identify its problems and how to fix them.

A while back, I worked on a similar problem involving inputting two time domain signals and classifying it (input: [# samples]-by-2, output: probabilities corresponding to each class). Now you will input the two time domain signals and you need to think carefully about how your output is formatted because sometimes you might have one activity, other times you might have five activities, or sometimes you might have no activity.

For example, say you input two signals ($x_1[n], x_2[n]$) of length 1024 and say you had activity #1 during samples 0-255, activity #2 from 256-511, and all other time is activity #0. Then you could format your output to be a signal of length 1024 (equal to length of input signals), to be:
$y[n]=
 \begin{cases}
      1, & 0 \leq n \leq 255 \\
      2, & 256 \leq n \leq 511 \\
      0, & n \geq 512
   \end{cases}
$. I bet this requires a bit of preprocessing of your dataset but this is not uncommon at all when doing these sorts of things.


Difference between repeated poles and distinct poles?

Repeated poles simply means there are more than one pole at the same location. If a pole is not repeated then it is a Distinct Pole. Consider the simple case of a cascade of two integrators (in s) which is similar to the cascade of two accumulators (in z).

The Laplace transform of an integrator
Is $\frac{1}{s}$, which is one pole at the origin. Two integrators in cascade would be $\frac{1}{s^2}$ so has repeated poles at the origin. The same with two accumulators in cascade would be $\frac{1}{(1-z^{-1})^2}$, which has repeated poles at z=1.

If there were two poles at z=1 and a third pole at z= 0.5, then such a system would have two distinct poles.


QPSK demodulator on GNURadio continuously loses phase lock

From what I understand, the modulation is a QPSK at 72k symbol.
So the main lobe is 144 kHz wide.
First point: Shaping filter clarification
Do you know if there is a shaping filter used by the transmitter?
Usually it is an RRC filter with a roll off of a given value.
If there is one in the transmiter, your RRC filter in the receiver should match roll off value
Second point: The low pass filter
The main lobe is 144kHz wide (because of QPSK-72ksps). The low pass filter you use is 90 kHz (2 x 45 kHz, for negative + postive spectrum side) wide.
With a 0.5 alpha roll off used at the transmitter, its bandwith would be 72 x (1+0.5) = 108 kHz <=> 54 x 2. With a 45 kHz, you cut before a possible RRC inside the chain (45 vs 54). Expand the low pass filter bandwith so that it does not interfere with the RRC filter in the receiver of Symbol block.
Try 60 instead of 45 Khz
Third point: AGC

AGC placement:
Is the AGC placement intended?
What happens if the signal is too low at the input of the first filter?
Do you have an external AGC before entering GNU radio?

AGC rate:
I have a small experience with GNU radio, so don't know what is behind the rate or other AGC blocks...
The important thing is that the AGC does not track on the varation of the symbols but on the average signal level.
Try maybe to decrease that rate so the reponse is slower.

There are probably other points to tune, these are main I see. You could earn several fractions of dB in those first places.
If there is an improve, you should see it on the constellation directly (if not too noisy)
The best approach is to first make the design and verify it performances in lab conditions and not directly on the satellite.


How do I know if my EEG signal need denoising?

It primarily relies on,

what are you looking for in the EEG: define the frequency band that
relates with the phenomena you are interested to investigate and
filter out the rest (e.g. for sleep events bandpass ~= 0.8-30Hz).
power line of your area: notch filter at 50Hz or 60Hz.
Study design: if you are investigating events ensemble averages will remove random noise processes and highlight potential EEG events of interest.

Moreover, you can apply baseline removal to de-noise your signal from saturation trends and other systemical artefacts. Other experimental methods include ICA, which you will not be able to apply robustly since you only have two channels.


Are location in BRIEF feature descriptor reused?

When using a randomized pattern in BRIEF, this means that you computed random positions inside the patch once in an offline procedure, then used these random locations every time you computed the descriptors. This makes sense, as it means that when comparing descriptors you will actually compare the same locations, it's simply that the sampling pattern was not obtained by some regular or geometric scheme (which is probably what makes it a bit more generic than hand-crafted patterns).
For example, here is the sampling pattern used by OpenCV.


You only get one row of three dots because the image is constant in the vertical direction.  That means that vertically speaking, the image only has 0 Hz content.  Thus, you get the three dots in row 0 (the 0 Hz row), and nothing in any of the others.  If you made one of the image's rows different than the others then you would see differences in the FFT's rows as well.


What is the meaning of channel capacity?

1) Is there a connection between the modulation kind and the channel capacity?

The capacity of a channel indicates the upper limit of how many bits can be transmitted per second over the channel with no errors (okay, technically it is "arbitrarily low number of errors", but it's basically the same thing).  We do various things to try to get as close to that upper limit as we can, such as modulation types, error correcting codes, etc., but none of them affect the channel capacity itself.

The equation for channel capacity is:

$$
C=B*\log_2(1+\frac{S}{N})
$$
where $B$ is the channel bandwidth in Hz, and $\frac{S}{N}$ is the signal to noise ratio (which is unitless).  Given this, the channel capacity only depends on its bandwidth and the received signal-to-noise ratio.

Going to higher order modulation types (e.g. from QPSK to 16-QAM, from 16-QAM to 256-QAM, etc.) is one way to take advantage of the $\frac{S}{N}$ portion of the equation.  If you have a very high signal-to-noise ratio it is wasteful to transmit QPSK, because the constellation points are very far apart from each other.  By going to 16-QAM the points are closer together, but they are still far enough apart to avoid most errors while transmitting twice as many bits.  If you still have "excessive" signal-to-noise you could go to a higher order QAM signal yet, allowing more bits per symbol to be transmitted.

  Why and when do we do $\log_2(M)$ where $M$ is the number of modulation symbols? For example, if the modulation is 16 QAM, then what does $\log_2(16)$ tell us?

$M$ is not the number of modulation symbols, it is the number of constellation points.  What $\log_2(M)$ tells us is the number of bits that each symbol represents.  There are 16 constellation points in 16-QAM, so we could number them from 0 to 15.  How many bits does it take to represent 0 to 15?  4 bits, or $\log_2(16)$.

  3) In equalization, when pilot symbols are send, we get better estimation quality. This is known as non-blind equalization. But in many research articles often it is mentioned that non-blind technique wastes a lot of bandwidth due to periodically sending the training symbols. What is the meaning of bandwidth wastage and how does it occur?

What they mean is that some of the bits that could have been used to transmit information are used to transmit pilot sequences instead.  In LTE, for instance, in some OFDM symbols every third subcarrier is a pilot instead of data.  If each subcarrier is using QPSK modulation, then each pilot subcarrier means that there are two fewer data bits.

  Does capacity decreases with using training/ pilot symbols?

No, not necessarily.  Yes, the number of information bits is decreased by the pilot symbols as described above, but remember that the channel capacity refers to how many error-free bits can be transmitted over the channel.  The pilots are essential for error-free communication.

A similar example is error correction codes (e.g. Turbo codes).  Error correction codes use parity bits or symbols to detect and correct errors.  We could replace those parity bits with information bits, but that would actually reduce the number of error-free bits that we sent through, so using the error correction code is a net gain.


Prove that time-reverse yields inverse filter for deconvolution?

This is not a general result. It specifically holds, iff the Fourier transform of the signal is of unit magnitude. In that case you can write the signal in frequency domain as

$$S(\omega)=\exp(i \phi(\omega))$$

for a real function $\phi$. The time reversed signal (I'll ignore the constant time offset, it doesn't change anything relevant) is then the conjugate

$$\bar{S}(\omega)=\exp(-i \phi(\omega))$$

Convolving the two by multiplying the frequency domain representation gives

$$\exp(i\phi(\omega))\times \exp(-i\phi(\omega)) = 1$$

Which is the response of the identity system and demonstrates that the signal and its time reversed are the responses of inverses.

This result is only exact in this special case. Sweeps can belong to this class of signal, for example if they have been generated as an allpass filter response. You can also do that by constructing the sweep entirely in frequency domain.

A sweep that is synthesised in time domain using a modulated sinusoid is typically not of that form, but approximates it well. Especially near the endpoints of the frequency range you will get significant deviation from the unit response.


How to convolve $u(-t)$ with other signals?

In general, convolution of a function $f(t)$ with $u(-t)$ results in a function $g(t)$ defined by the following integral:

$$g(t)=\int_{-\infty}^{\infty}f(\tau)u(\tau-t)d\tau=\int_{t}^{\infty}f(\tau)d\tau\tag{1}$$

In the given example, $f(t)=u(t+1)$, so you get

$$g(t)=\int_{t}^{\infty}u(\tau+1)d\tau=\int_{t+1}^{\infty}u(\tau)d\tau\tag{2}$$

I leave it up to you to figure out what the result of that integral is.


SVD Based Noise Filter on the Samples in Toeplitz Form

The matrix M is called a Toeplitz matrix and they have a certain type of decomposition (but X should be the ACF and not the signal itself), so the resulting U and V would be a Vandermonde matrix which are very similar to Discrete Fourier transform.
So the whole process is similar to find spectrum of your signal and putting zero in bins which their power are below a threshold.


When you filter a finite length signal, the filtered output will have a larger length, either you implement the filter in time domain or frequency domain. Your desired continuity comes from this excess part which must be added to the result of next segment or window.

For example when you take N samples of a signal and filter it with a filter of length M, the result has a length of N+M-1, where the M-1 excessive samples must be added to the next segment. To do so you have to find your filter response in time domain which has a length of M samples, then zero-pad both your signal and filter to N+M-1 samples then take N+M-1 point FT of both and multiply them.


Is this graph of Wigner-Ville distribution OK?

The graph is over amplified. If the image in blue is divided by factor of 20, the transients vanish and only a line of the detected frequency is visible.


How do I test stability of a MIMO system?

Let's try to make your system into a state space representation with state:
$$
\mathbf{x}(t) = \left [
\begin{array}{c}
x_1(t)\\
x_2(t)\\
x_1(t-1)\\
x_2(t-1)\\
\end{array}
\right ]
$$
so that
$$
\mathbf{x}(t) = \mathbf{A} \mathbf{x}(t-1) + \mathbf{y}(t)
$$
where
$$
\mathbf{A} =
\left [ \begin{array}{cccc}
a_{11} & a_{12} & a_{13} & a_{14}\\
a_{21} & a_{22} & a_{23} & a_{24}\\
1   & 0 & 0 & 0\\
0   & 1 & 0 & 0\\
\end{array}
\right ]
$$
and
$$
\mathbf{y}(t) = \left [
\begin{array}{c}
y_1(t)\\
y_2(t)\\
0\\
0
\end{array}
\right ]
$$
or you could make $\mathbf{y}$ $2 \times 1$ and introduce a $\mathbf{B}$ matrix.

The system may be unstable if the eigenvalues of $\mathbf{A}$ are outside the unit circle.


What is the purpose of wrapping the negative times of a response function in discrete convolution?

It's about zero-phase filtering.
First, the continuous-time response function $r(t)$ is an even-symmetric, noncausal filter so that its Fourier transform will be a zero-phase function; i.e.,
$R(\omega)$ is a real-valued function of frequency $\omega$.
Note that there's almost a mistake in the figure 13.1.1; the output should begin before $t=0$, even for a zero-phase filter. Because this is a non-causal filter, it will anticipate future inputs before they happen, and put out the response before $t=0$.
I think, however, they wanted to justify this condition, by shifting the input $s(t)$ sufficiently to the right (or assuming that $s(t)=0$ until $t=t_1$) so that the output response before $t=0$ will be zero anyway, and thus eventually we are interested at the output of the filter beginning from $t=0$.
Now, sample the signals $s(t)$ and $r(t)$ (assuming they are bandlimited) and obtain a sampled version of the continuous-time output by performing a discrete-time linear convolution, or equivalently use a DFT based efficient implementation of linear convolution by circular convolution of proper length.
If you choose the first approach, and perform direct discrete-time linear convolution using the sequences $s[n]$ and $r[n]$, then you won't need to shift anything. But remember that a practical implementation will still need some shifts to represent the noncausal impulse response $r[n]$ in the proper index frame.
But if you choose the second approach; using DFT to implement a linear filtering, you have two choices for the filter response $r[n]$ to be registered; (because in the DFT domain, there's no negative time-index $n$, and you are restricted to the range $0 \leq n < N$) first you can linearly shift $r[n]$ sufficiently right to make it a causal filter $r_c[n]$, and follow the procedure as usual; with this choice, the effective discrete-time filter loses its zero-phase (and zero group delay) property and your effective output begins at the group delay position index.
If, instead, you wrap the negative-index portion of $r[n]$ to its far right edge, in the DFT frame of $N$ samples, where $N$ is chosen properly to aviod aliasing in the time-domain to ensure that computed circular convolution will match the desired linear one, then you will maintain the zero-phase property of $r[n]$ while performing the DFT based convolution; your first output at $n=0$ will be what you would get at $t=0$ in the continuous-time case.


What is the cut-off frequency of a moving average filter?

The moving average filter (sometimes known colloquially as a boxcar filter) has a rectangular impulse response:

$$
h[n] = \frac{1}{N}\sum_{k=0}^{N-1} \delta[n-k]
$$

Or, stated differently:

$$
h[n] = \begin{cases}
\frac{1}{N}, && 0 \le n < N \\
0, && \text{otherwise}
\end{cases}
$$

Remembering that a discrete-time system's frequency response is equal to the discrete-time Fourier transform of its impulse response, we can calculate it as follows:

$$
\begin{align}
H(\omega) &= \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} \\
&= \frac{1}{N}\sum_{n=0}^{N-1} e^{-j\omega n}
\end{align}
$$

To simplify this, we can use the known formula for the sum of the first $N$ terms of a geometric series:

$$
\sum_{n=0}^{N-1} e^{-j\omega n} = \frac{1-e^{-j \omega N}}{1 - e^{-j\omega}}
$$

What we're most interested in for your case is the magnitude response of the filter, $|H(\omega)|$. Using a couple simple manipulations, we can get that in an easier-to-comprehend form:

$$
\begin{align}
H(\omega) &= \frac{1}{N}\sum_{n=0}^{N-1} e^{-j\omega n} \\
&= \frac{1}{N} \frac{1-e^{-j \omega N}}{1 - e^{-j\omega}} \\
&= \frac{1}{N} \frac{e^{-j \omega N/2}}{e^{-j \omega/2}} \frac{e^{j\omega N/2} - e^{-j\omega N/2}}{e^{j\omega /2} - e^{-j\omega /2}}
\end{align}
$$

This may not look any easier to understand. However, due to Euler's identity, recall that:

$$
\sin(\omega) = \frac{e^{j\omega} - e^{-j\omega}}{j2}
$$

Therefore, we can write the above as:

$$
\begin{align}
H(\omega) &= \frac{1}{N} \frac{e^{-j \omega N/2}}{e^{-j \omega/2}} \frac{j2 \sin\left(\frac{\omega N}{2}\right)}{j2 \sin\left(\frac{\omega}{2}\right)} \\
&= \frac{1}{N} \frac{e^{-j \omega N/2}}{e^{-j \omega/2}} \frac{\sin\left(\frac{\omega N}{2}\right)}{\sin\left(\frac{\omega}{2}\right)}
\end{align}
$$

As I stated before, what you're really concerned about is the magnitude of the frequency response. So, we can take the magnitude of the above to simplify it further:

$$
|H(\omega)| = \frac{1}{N} \left|\frac{\sin\left(\frac{\omega N}{2}\right)}{\sin\left(\frac{\omega}{2}\right)}\right|
$$

Note: We are able to drop the exponential terms out because they don't affect the magnitude of the result; $|e^{j\omega}| = 1$ for all values of $\omega$. Since $|xy| = |x||y|$ for any two finite complex numbers $x$ and $y$, we can conclude that the presence of the exponential terms don't affect the overall magnitude response (instead, they affect the system's phase response).

The resulting function inside the magnitude brackets is a form of a Dirichlet kernel. It is sometimes called a periodic sinc function, because it resembles the sinc function somewhat in appearance, but is periodic instead.

Anyway, since the definition of cutoff frequency is somewhat underspecified (-3 dB point? -6 dB point? first sidelobe null?), you can use the above equation to solve for whatever you need. Specifically, you can do the following:

Set $|H(\omega)|$ to the value corresponding to the filter response that you want at the cutoff frequency.
Set $\omega$ equal to the cutoff frequency. To map a continuous-time frequency to the discrete-time domain, remember that $\omega = 2\pi \frac{f}{f_s}$, where $f_s$ is your sample rate.
Find the value of $N$ that gives you the best agreement between the left and right hand sides of the equation. That should be the length of your moving average.


With regards to multirate processing, what are the benefits of having a 'slower' sampled signal?

In addition to what hotpaw2 said, it's important to understand that sample rate reduction is the corner stone of actual applicability in many systems.

Software Defined Radio frontends nowadays produce hundreds of complex-valued Megasamples per second – if you don't reduce that amount of data, no commonly available PC can deal with that.

Also, many devices that you meet nowadays actually have software-programmable DSP embedded very deeply. The average mobile phone baseband IC, i.e. what handles the signal, is digital and very broadband, to "get" the whole bandwidth that e.g. LTE demands they can capture; however, at a single time, only a fraction of that bandwidth might be relevant. So isolating/extracting that bandwidth is crucial to functionality.

Not reducing the sample rate in that situation would not only be a pure waste of power – it would render the following steps, such as equalization and synchronization, totally impractical in the telephone form factor.

A movie camera might take up to let's say 16 Million pixels – which are, in the end, nothing but samples of a twodimensional signal! If the producing camera didn't remove that sample rate (i.e. the picture's resolution), there's no way that the resulting movie data would fit through the terrestrial TV channel; no way that a cheap TV might decode it, and no way you could by such a movie on a cheap optical disc at current capacities.

So, this, by far, is not a "historical" problem. It is a very current problem – not being able to handle the full sampling rate that ADC/DAC hardware offers is the core technological problem of pushing communication further.

Now, 10 kHz really isn't something you'd normally worry about – that's something that most algorithms could deal with on the cheapest microcontroller. But think more in the terms of a SDR frontend (e.g. this one) that has 200 MHz of sampling rate per single complex ADC; that means that you can't even get the  two channels of that device into whatever computer you attach at the native sampling rate, because $2\cdot200 \cdot 10^6 \frac {\text S}{\text s} \cdot \frac{16 \text{ b [I]} +16 \text{ b [I]} }{\text S} = 12.8 \frac{\text{Gb}}{\text s} > 10 \frac{\text{Gb}}{\text s}  $, which is the datacenter-grade 10GbitEthernet interface's data rate that usually totally flattens out even the most powerful quad-core PC CPUs, without them even doing any DSP, just shuffling the data into RAM. Now, you realize that for almost all applications, it's very desirable that the device internally reduces the sampling rate of the observed signal as far as possible to correctly represent that signal, for computers to even be able to cope with the raw amount of data. The device does that, by applying cascaded multirate filters internally. The ability to, nevertheless, observe the whole bandwidth at once and "select" any part of these 200 MHz at random is very crucial to designing things like wideband transceivers or frequency-hopping standards; so the ability to use the minimum viable sampling rate is very central to progress in research, application development and implementation, as soon as you're not considering 10 kHz systems any more.


Any Open Source Fast Wavelet transform libraries?

you can have a look at the LTFAT's wavelet module
http://ltfat.sourceforge.net/doc/wavelets/index.php

it runs in Matlab/Octave with backend written in C. It has fairly large database of wavelet filters and new ones can be added easily.

What exactly do you mean by

  2) Ability to run in parallel - VERY IMPORTANT

Should the computation itself be somewhat parallelized or that it should be possible to run several batches of FWT calculations? I guess that running several instances of Octave should do the trick in the latter case.


Equalization : using the spectral or the temporal signal?

If you are trying to make a proper filter, you would want to use FFT convolution, which is like OLA, but not the same.  OLA is more of a synthesis technique for reducing noise between frames caused by incoherent phase mangling from directly manipulating frequency data.
You also need to resolve the phase of the bins for your filter.  You are only defining the magnitudes and phase at those discrete frequencies, not the frequencies in between, which can be radically different then what you have in your head.  The easiest thing to do is to use linear phase to center the impulse of the filter in the window.  This can be done by multiplying each bin magnitude by e^jpik.
If you are trying to realize an analog prototype, this is a bad way to do it.  Your frequency bins are linearly spaced, but analog filters are logarithmic.  The effect will be that your lower bands have a lower “Q” and vice versa.  You can still do it, but it may not do what you wanted.


The fundamental property of averaging filters (in swear words, a linear system) is that the output of a sine is a sine of the same frequency, albeit of zero amplitude. So:

Which non-trivial moving average weights would you need to find to get a filter that returns the same output as its input?

An averaging in phase with the sine: take a filter full of zeroes, except $K$ values, pairwise distant by a multiple of $40$, with coefficient values $1/K$. This way, you will average $K$ sine value of the same amplitude, yielding a value with the same amplitude.

Hence, you have a bunch of non trivial filters. However, they are not very interesting in practice, because they shall be tuned to the period of a sine sampled with particular sampling.


You're right, the given $x[n]$ is clearly periodic. You can show this by simply checking if

$$x[n]\stackrel{?}{=}x[n+N]\tag{1}$$

is satisfied for some positive integer $N$.

For the given $x[n]$ you get

$$\begin{align}(-1)^{n^2}&\stackrel{?}{=}(-1)^{(n+N)^2}\\&=(-1)^{n^2}(-1)^{N(N+2n)}\tag{2}\end{align}$$

From $(2)$ it follows that for $(1)$ (and $(2)$) to hold we require

$$(-1)^{N(N+2n)}=1\tag{3}$$

i.e., $N(N+2n)$ must be even (for arbitrary integer $n$). This is certainly the case for any even $N$, and the smallest positive even $N$ is $2$, hence the (fundamental) period of the given $x[n]$ is indeed $2$.