Relation between Fourier phase and Analytic Signal Phase

Discussing the relationship between the Fourier transform and the Hilbert transform, let us consider an example function. It can be any square integrable function; for readability of the text, we select a function with both transforms expressible via simple analytical functions.
Let a signal s be a sinc function of time, s(t) = sinc(t). The Fourier transform of the sinc function is a rect function
$$
\mathscr{F}\{s(t)\}(ω)={{1\over{\sqrt{2π}}}·rect\left({ω\over{2π}}\right)}
$$
which is a function of frequency: OP's FT_s is a function of frequency, or, for the discrete case, an array with indices representing frequency bins.
The Hilbert transform of the sinc function
$$
\mathscr{H}\{s(t)\}(t)={{1-cos(t)}\over{t}}
$$
is the function of time, because the Hilbert transform is not a transform of the kind that maps time domain functions to frequency domain functions -- as the Fourier transform does -- but "a specific linear operator that takes a function, u(t) of a real variable and produces another function of a real variable H(u)(t)" (citation from the Wikipedia article on the Hilbert transform).
Comparing these formulas, one sees that the question "if there is a relation between the phase of the Fourier Transform and the phase of the Analytic signal?" is rather pointless.
FT_s, as defined by OP, is an array of values against frequency bins, and atan(Imag(FT_s),Real(FT_s)) is also an array indexed by the same bin enumerables, corresponding to a continuous-time frequency. Therefore, atan(Imag(FT_s),Real(FT_s)) is an array of phases of the signal frequency components.
A_s is a function of time, and atan(Imag(A_s),Real(A_s)) is also the function of time. It is the instantaneous phase of an analytic signal, and it is not the phase of the Fourier transform frequency component, not in any way.
The OP-defined two continuous-variable "phase functions" and corresponding discrete "phase arrays", phaseFourier and phaseAnalytic, are defined on different domains: phaseFourier is defined on a frequency domain, phaseAnalyticis defined on a time domain. Still, you may want to see the graphs of these function for a sinc(t) signal.
For s(t) = sinc(t), the Fourier transform (OP's FT_s) values is a real constant up to a frequency of 1/(2π) and a zero for higher frequencies. Therefore, the array of phases of frequency components (OP's phaseFourier) is zero up to 1/(2π) and undefined in the region where the Fourier transform amplitude is zero. For sinc(t), one can safely assign all phases of FT components to zero in the entire frequency domain. The FT_s graph is trivial, and shown in every textbook and course notes on signal processing.
The instantaneous phase (OP's phaseAnalytic) graph may be shown less frequently, so here you are:

Summing up, one may agree that, talking about the relationship between the Fourier transform and the Hilbert transform, more appropriate is to indicate that the Hilbert transform is a multiplier operator. The multiplier of H is
$$
σ_{\mathscr{H}}(ω)=-i·sgn(ω)
$$
(see https://en.wikipedia.org/wiki/Hilbert_transform#Relationship_with_the_Fourier_transform).


The only thing you need to know is that an LTI system's response to a sinusoidal input signal is a sinusoidal output signal with the same frequency as the input signal, with its amplitude and phase modified by the magnitude and phase of the system's transfer function evaluated at the input frequency, or, in math:

$$x[n]=A\cos(\omega_0n)\Longrightarrow y[n]=A|H(\omega_0)|\cos(\omega_0n+\phi(\omega_0))\tag{1}$$

where $x[n]$ and $y[n]$ are the input and output signals, respectively, and $H(\omega)=|H(\omega)|e^{j\phi(\omega)}$ is the system's frequency response.

So if you got $H(\omega)$, you just need to evaluate it at the given input frequency $\omega_0=\pi/2$ and apply $(1)$. Splitting up the cosine into two complex exponentials, as suggested by the given solution, is not necessary (and not even easier I'd say).


Perform a 2 dimensional cross correlation against a noise free replica of the signal you expect to receive with delay offset on one axis and frequency offset on the other, and find your maximum correlation. Normalize that to get $\rho$ and from $\rho$ you can get SNR. Do this after channel equalization, otherwise you will get the SNR degradation due to ISI as well.

For further details, see:

Noise detection

and

How can I find SNR, PEAQ, and ODG values by comparing two audios?

If the SNR is high and the noise is well behaved (thermal noise floor, no band selective interference present) then you can make a reasonable estimate directly from a spectrum plot such as the example I show below, where the SNR in this case was 40 dB as is clear from the distance between the signal and the noise. This is signal power to the total noise power that is in the signal bandwidth.


There are two versions of optical flow(OF): Feature based (sparse) or dense. In the dense version OF is applied to all the image pixels, while in the sparse one, only certain characteristic feature points are tracked. However, both approaches depend on the tracking of pixel quantities. This is fundamentally different than tracking the whole patch, because in return one obtains a full set of pixel-level correspondences. This is not just a vague statement, but rather carries the idea that the correspondence estimation is not constrained. Remember in the motion case, one would assume that the motion is "rigid" or "articulated" and imposes this prior into the tracking framework. If optical flow is used to estimate the dense trajectories, than this constraint is not assumed and one could as well track the deformable bodies, regardless of the deformation model.

Optical flow provides you more freedom and information about the tracked scene. Yet, my experience is that the motion tracking methods are more robust and reliable. But keep in mind that, by applying human models or temporal models on top of such tracked trajectories (from OF) researchers are capable of developing robust articulated / deformable / model-based tracking algorithms.


FFT of image data: "mirroring" to avoid boundary effects

The FFT of the original image is correct. The artifacts you are seeing are typical for the DFT, because the DFT basis functions have difficulty representing non-periodic signals. Even though the DFT is of finite length, it actually represents a periodic signal that extends to negative and positive infinity. The DFT's basis functions are all sinusoids with periods that are integer divisors of the DFT size, so they all start and end at the same value, and it is hard for it to fit signals that are not periodic in this way. So, imagine tiling your image: there is indeed an abrupt discontinuity at the edges. The edges do not match up well, and that is the reason for the boundary effects.

In audio signal processing, where the FFT is used often, a very common approach is to use a window function, as you noted. This tapers the edges off towards 0, and reduces this effect.

In image processing, the DCT is usually used rather than the DFT, because it imposes different symmetry constraints that give better results. In other words, it's used to prevent exactly the issue you're seeing. In contrast to the DFT's basis functions (remember how they all start and end at the same value), the DCT's basis functions are integer divisors or half of integer divisors, so many of them start and end at different values. As a result, the implied boundary conditions are different: the signal is assumed to be symmetric about its edges, which is actually somewhat similar to the "mirroring" you were experimenting with. Here is another post by me with a little more information on the DCT: https://dsp.stackexchange.com/a/362/392

So, as a short summary, if you are absolutely set on using the FFT, you might want to try windowing. However, the best choice is probably to use the DCT, which implies boundary conditions similar to your "mirroring" idea, and as a result, handles images better.


,,A wavelet'' actually is two functions that are orthogonal to each other. That way, a signal can be split into two parts. One part is the wavelet proper, while the other is usually called ,,scaling function'' (sometimes ,,mother wavelet'' as well). When you have just one of both, you can not perform a meaningful wavelet transform.


How to recognize for small video sequence with a large sequence?

The most intuitive approach here would simply be taking every pixel's color values, and finding peaks in the individual pixels' reference video/short video cross correlation function. That's basically a W×H-dimensional crosscorrelation function; you could speed that up by doing fast convolution.

You could algorithmically reduce the effort by first reducing the dimensionality of your problem, for example by reducing each video frame to a global brightness value and a mean motion vector of e.g. the frame center region – then you'd only have to calculate a threedimensional crosscorellation function.


In your images, it seems regions are defined more by their edge boundaries rather than their colour.

Since the IFT is using a colour difference metric then it doesn't see much difference between some foreground and background super pixels.

Instead of colour difference I would try a different metric based on boundary edge strength.   For example, the weight of the edge between two nodes is set to the average edge strength along their shared super pixel boundary. You might need to define a minimum boundary length though, such as 20 pixels, otherwise noise can effect this.

I've used this measurement before to combine segmentation regions and it works well. Interestingly, I was getting the segmentations using the IFT on a per pixel basis. This method you are applying looks interesting!

You could also combine the two measurements.


How to detect time of arousal with excitation signal in speech signal?

You are looking for the instants of glottal closure, or GCIs.  There are many freely Matlab algorithms to do such.  Check out, for example, the DYPSA algorithm.  A paper explaining how these algorithms generally work can be found here.


Finding the Width of a Point Spread Function (PSF)

One simple way is to identify zeros in the Frequency Domain.
This was one of the legitimate method for Blind Deconvolution few decades ago.

Draw the connection between the zeros in frequency domain of the Convolution Kernel to the parameter $ w $.
Learn the zeros in frequency domain of the given image.
Find the $ w $ which explains most of the zeros of the image.

It's pretty simple and for a parameterized family of PSF's can be quite efficient.


How to perform dot annotation for crowd counting in a image file using Matlab?

The Viola Jones algorithm does the face detection for you using the Cascade Object Detector.

%Read the image
x= imread('C:\~\crowd.jpeg');
x=rgb2gray(x);

%Face detection
faceDetector = vision.CascadeObjectDetector;
bboxes = step(faceDetector,x);

The first two columns of bboxes gives the (x,y) coordinates of the faces detected

If you want to label the people in the picture:

x1 = bboxes(:,1) %x co-ordinates
x2 = bboxes(:,2) %y co-ordinates
figure;
imshow(x);
for i=1:size(x1,1)
    plot(x1(i),x2(i),'*');
end

Saving the data in a .mat format:

save('dot.mat','x1','x2');

Accessing the data:

load('dot.mat');

To save data points from a given image, you can use ginput method. First read the image and then use [x,y] = ginput(n) where n is the number of points you need to save.


Is there a closed form expression for noise other than white noise?

This is usually done with filtering of white noise.

White noise generator + LP filter give you a red noise generator.
White noise generator + HP filter give you a blue noise generator

Design a filter with required frequency response and apply it to a white noise generator.


Linear Constant Coefficient Differential Equations: Zero-Input and Zero-State responses

The two approaches should return the same solution. They are just two different ways to get to the same place.

In the ZSR/ZIR method, you are solving two different IVPs - they have the same expression on the left-side of the equality, but different initial conditions and different right-hand sides (one has a $0$, the another one has some function of $t$). In the more mathematical non-engineer like approach (homogeneous/particular), you solve just one IVP with both things: the right-hand side of the equation not set to $0$ and the original initial conditions. Thus one can intuitively see the second approach as some kind of superposition of the two IVPs proposed in the first one.

Just in case you don't believe this yet, let me just write down an example (without doing it in detail as that's not the point of this post) to show that both methods work. Let's solve the following differential equation:

$$\left\{
 \begin{array}{ll}
  y''(t)+4y'(t)+3y(t)=e^{-2t} \\
  y(0)=1 \ , \ y'(0)=1
 \end{array}
\right.$$

The mathematical approach ($y_p+y_h$)

I propose:

$$y_p(t)=-e^{-2t}$$

Check that it satisfies the differential equation and the initial conditions:

$$y_p''(t)+4y_p'(t)+3y_p(t)=-4e^{-2t}+8e^{-2t}-3e^{-2t}=e^{-2t} \ \implies\mathrm{OK}$$
$$y_p(0)=-e^0=1\ \implies\mathrm{OK}$$
$$y_p'(0)=-2e^0=1\ \implies\mathrm{OK}$$

Therefore we conclude that $y_p(t)$ is a valid particular solution. Then we should solve the homogeneous equation. I'm not going to do it step by step, but you can check thay you get something of the form $$y_h(t)=Ae^{-t}+Be^{-3t}$$

If we impose the initial conditions to $y_p(t)+y_h(t)$, the constants turn out to be $A=5/2$ and $B=-1/2$. So the solution can be written as:

$$y(t)=y_p(t)+y_h(t)=\frac52 e^{-t}-e^{-2t}-\frac12 e^{-3t}$$

The engineer approach ($y_{zs}+y_{zi}$)

As stated at the beginning, this approach requires us to solve two different IVPs. The first one is:

$$\left\{
 \begin{array}{ll}
  y_{zi}''(t)+4y_{zi}'(t)+3y_{zi}(t)=0 \\
  y_{zi}(0)=1 \ , \ y_{zi}'(0)=1
 \end{array}
\right.$$

The solution (the generic one, with the constants $A$ and $B$) to this problem is the same as in the previous case, but with different constants (now we only have the homogeneous solution term). In this case, $A=2$ and $B=-1$. Therefore:

$$y_{zi}(t)=2e^{-t}-e^{-3t}$$

Now we have the other IVP:

$$\left\{
 \begin{array}{ll}
  y_{zs}''(t)+4y_{zs}'(t)+3y_{zs}(t)=e^{-2t} \\
  y_{zs}(0)=0 \ , \ y_{zs}'(0)=0
 \end{array}
\right.$$

Note that his is exactly the same IVP as the one stated at the beginning, but with different initial conditions. So we can use the homogeneous/particular approach to solve this IVP as well, and find the constants $A$ and $B$ that suit the initial conditions in this case. We find that

$$y_{zs}(t)=\frac12 e^{-t}-2e^{-2t}+\frac12 e^{-3t}$$

If we sum the two solutions:

$$y(t)=y_{zs}(t)+y_{zi}(t)=\frac52 e^{-t}-e^{-2t}-\frac12 e^{-3t}$$

We got the same result.

Conclusions

As can be seen from the example, both methods arrive to the same solution because they are just two different ways of decomposing the problem into two easier ones, but they are equivalent.

The approach using ZSR and ZIR is more practical in real life. Take for example a circuit. We usually know the initial conditions with no input. When we turn the input on (i.e. we make a sudden change in the circuit), voltages and currents could change abruptly. The ZIR uses these intitial conditions with no input, and the ZSR uses the input when we know all initial conditions are zero (for example, we discharge all capacitors before turning the circuit on). Performing the mathematical approach would not be simple, as we would have to know the initial conditions just when the input appeared, and that can be rather hard to measure or observe. With ZSR/ZIR, we can measure the two things independently and then add them up to get the total solution.

Regarding the steady state question, note that when you solve the differential equation you get the signal $y(t)$. Namely, you can evaluate the output at any instant $t$ you want. Thus, you don't have just the transient response or the steady state one: you have both. You can observe whatever you want as you have the signal as a function of time. The steady state will be reached at $t\to\infty$ (in real life, after some time, when all components in the circuit are acting as the input had been turned on for a long time).


Basically you want to implement a filter with batch processing, right?

Is this an FIR filter or an IIR filter?

What is the sampling rate? Number of taps? Acceptable latency between output and input as batch processing adds some processing latency.


Find the autocorrelation function of signal $x(t) = u(t) - u(t-1)$

Don't make this more complicated than it really is. $x(t)$ is non-zero in the interval $t\in[0,1]$, and $x(t+\tau)$ is non-zero in the interval $t\in[-\tau,1-\tau]$. The integrand is non-zero only if the two functions overlap. There is no overlap for $1-\tau<0$ and $-\tau>1$, i.e., for $|\tau|>1$. So for $|\tau|>1$ the autocorrelation is zero.
For $-1<\tau<1$, we have, according to above considerations, the following integral:
$$\mathcal{R}_{xx}(\tau)=\int_{\max\{0,-\tau\}}^{\min\{1,1-\tau\}}dt=\begin{cases}\displaystyle\int_{0}^{1-\tau}dt,&0<\tau<1\\\displaystyle\int_{-\tau}^{1}dt,&-1<\tau<0\end{cases}$$


Is there a faster algorithm about this type of FFT

Alas, a multiplication by a rectangular window in the time domain is a convolution in the frequency domain.  Thus (assuming no other constraints on or knowledge about A or B), the DFT of A (or B) would have to be the convolution of the DFT of C with a twisty Sinc function, which is computationally far more expensive than 2 FFTs (except by fast convolution, which would require the same number of FFTs as your "ordinary" method).


Larger image after wavelet reconstruction

A generic method consists in adapting the wavelet decomposition on the border, ie the wavelet varies on shorter frames. It can be called "wavelet on the interval" or be addressed with lifting schemes. But this requires caution and complicated computations sometimes.

A more lazy option (if I remember well)  is to use a pair of symmetric/anti-symmetric $M$-band wavelets, over $J$ levels. If you expand the image to a size that can be divided by $J^M$ (the most symmetrically as you can), and use  a symmetric extension that may reveal non-expansive if properly tuned. This may vary with implementations.


Does horizontal motion blur imply the shape of PSF is rectangular pulse?

Think a minute of what a motion blur means for the physical pixels: they are moved during the acquisition process.
Assuming left-to-right horizontal motion for simplicity, then what happens is that one pixel will measure the amount of light in its initial position, then the light in the position that was occupied before by its right neighbour, and so on until the frame acquisition ends.

Mathematically, the light measured by one pixel is then a mixture between the value in its original position and the values of its neighbours in their original positions.
Hence, up to a rescaling of the output value (that is unnecessary of you normalize the blur kernel), an horizontal motion blur can be modelled by an horizontal blur kernel.

Of course, this is grossly simplified (I did intentionally forgot about minimum exposure time for example and supposed instead that the measurement was instantaneous) that assumes that the scene is at a fixed distance from the camera and that all the image points are at the same distance, without independantly moving objects.


If my sampling rate is not high enough to capture all frequencies, can I accurately capture low frequencies?

If there are higher frequencies present in the signal you will indeed get all these! The problem is the frequency components will be aliased as additional signals also within your first Nyquist zone of 0-500 Hz (assuming a real signal, if complex then the full 0 to 1000 Hz spectrum in the DFT will be unique).

To uniquely observe the signals that are in the 0-500 Hz frequency range, you MUST have an anti-alias filter in the analog domain prior to the sampling of the signal at 1000 Hz. Ideally (but not achievable) you would use a "brick-wall" filter to pass the DC-500Hz band with no distortion or loss, and then completely reject everything above 500 Hz.  Welcome to the real world, where filters have finite impulse responses and we don't have forever to wait for long filter delays-- meaning we must have a transition band in the design of these filters. The amount of effort and cost you want to put into the anti-alias filter dictates how much of the DC-500 Hz band you can reliably observe, and really goes into a lot of consideration in the overall system design (trading design detail elsewhere). A typical number can be 35-40% (35% Meaning DC-350 Hz passband, 650 Hz + stop band.).

Now if your goal is to instead observe the spectrum that is present in the 1500 - 2500 Hz frequency band, you could instead design band pass anti-alias solutions to select them specifically prior to sampling with the same 1000 Hz sampling clock (this is called "under-sampling). In this case you have restricted what will alias into the DC-500 Hz band, and can therefore uniquely identify the actual frequency components from mapping the aliasing. Like above, if it is a real signal and you are restricted to the 1000 Hz sampling rate, then 350 MHz bandwidths would be reasonable for the bandpass filter, meaning you would not be able to observe with one DFT the full 1500-2500 Hz analog spectrum with a real sampling rate of 1000 Hz.

Undersampling is best demonstrated with the graphic below showing two analog signals, an 11 Hz sine wave and a 1 Hz sine wave both sampled with a 10 Hz sampling clock. Once sampled we will only have the red circles representing the signal. This shows that the input in this case could have been either the 11 Hz or the 1 Hz signal, and if no filtering was present to restrict one from being possible, we would have no way to distinguish which signal was present. If both are present then they could destructively or constructively interfere depending on their relative phase.


Compute absolute power spectrum in MATLAB

If you have the Signal Processing Toolbox, you can use pwelch. If you do >>help pwelch   There will be some other choices. Try a few. Only you have your data so only you can determine if your results are acceptable

If you don't have the Toolbox, using FFT isn't hard. It's a matter of the length of your data, window choice, overlap, and a few other parameters


Envelope of a signal via Convex hull / Edge detection

You can use the Hilbert transform to get the analytic signal.  The absolute value of the analytic signal is the envelope.  It's available in SciPy and many other software suites, though it's usually poorly named.


Difference equation and Impulse response of an IIR filter - deconvolution

So your transfer function is

$$H(z)=\frac{1}{1+az^{-L}}\tag{1}$$

with $a=0.25$ and $L=2205$. If $Y(z)$ and $X(z)$ are the $\mathcal{Z}$-transforms of the output and input signal of the filter, respectively, then you have

$$Y(z)=H(z)X(z)$$

which, using (1), is equivalent to

$$Y(z)(1+az^{-L})=X(z)\tag{2}$$

Transforming (2) back to the time domain results in the difference equation of the recursive filter:

$$y[n]+ay[n-L]=x[n]\tag{3}$$

You can get the corresponding impulse response either from the inverse $\mathcal{Z}$-transform of (1), or from setting $x[n]=\delta[n]$ in (3) and observing what happens to the output. In either case, the result should be

$$h[n]=\begin{cases}(-a)^{n/L},&& n=0,L,2L,3L,\ldots\\
0,&& \text{otherwise}\end{cases}$$


From an ADC perspective, it is just taking a sample of the voltage in time. I fail to see how a "misinterpretation" could be made since there is no "turning car wheel" to take pictures of at the wrong time. Do the harmonics alias in such a way that the wave shape is preserved?

You can reason this out yourself, in the time domain.  Consider a square wave with values -1 and 1; and sample it at exactly five times its period.  You'll get something like {+1, +1, -1, -1, -1, +1, +1 ...}.

If you just have to think about this in the frequency domain, you can, after a great deal of work, demonstrate to yourself that not only do the harmonics not alias in a way that preserve the shape, but in fact alias in a way that does not preserve the shape -- you end up with edges aligned with the ADC sample points.

  In my mind, adding a filter to the signal will modify the shape of the original signal, getting rid of those upper harmonics. Depending on the filter design, it could add overshoot, ripple, and/or rise/fall time changes. So wouldn't the best representation of the pulse be obtained by direct sampling?

That depends on what you're trying to do.  If your goal is to accurately capture the timing of the edges, to something less than the ADC sampling interval, then you need to round those edges out -- because a signal that goes {+1, +1, -1, -1} has less information about the location of the edge than one that goes, e.g. {+1, 0.75, -0.25, -1} or {+1, 0, -1, -1}.

And finally:

  I'm trying to determine whether or not an anti-alias filter is needed for sampling square waves.

That depends entirely on your problem at hand.  If you need every bit of information contained in the wave (and not just its timing), and if it truly has energy content out to infinity (which is physically impossible) then you need to sample infinitely fast.  If you know it's bandlimited, then you know the sampling rate you need to use.  If you know that, for example, the edges are super-sharp but what you really care about is their timing, then you can filter and acquire with an ADC, and infer the actual position of the edges from the way that the measured "square-ish" wave transitions.


Phase error correction for Fourier transform basis vectors

This can be expressed with vectors and matrices as $f'(x) = M\cdot f(\nu)$.  Here, the rows of $M$ correspond to an orthogonal set of sinusoidally varying basis vectors.

Woah, there, that's a big step you're doing. So, you have your band-limited signal; it has to be band-limited, otherwise the $x_\text{max}$ limit in $(1)$ would be illegal.

So, since it's limited in Fourier ($\nu$) domain, it has to have infinite support in time domain ($x$). (Imagine you take an infinitely long signal, and window it with a rectangular function: You need to convolve the Fourier transform of that spectrum with a scaled $\text{sinc}$ function, and that has infinite support)

Now, you're right to consider the Fourier transform as a linear operator $M$. What is very hard to do is justify that $M$ is actually a matrix with discrete rows of discrete frequencies. That integral there doesn't give you discrete, separate sinusoids, it gives you a continuum of frequencies!

Basically, the whole of the sampling theory is to prove that you can do that. Short idea: Thanks to Nyquist we know that we can represent any band-limited signal by sampling that signal at a sufficiently high rate, giving us, indeed, a discrete vector. The Time-Discrete Fourier transform of a signal collapses to a sampling rate-periodic spectrum, and we can show that this is also representable by a discrete vector, leading us to an understanding of the Discrete Fourier Transform (DFT), which actually is just a Matrix with complex numbers in it, of size $N\times N$, $N$ being the length of the time observation vector.

But: that discretization into an actual matrix over $\mathbb C$ requires you to first convert to discrete signals. Then, $M$ must be square, and is indeed (easily!) invertible. So, your problem lies in that somehow, you ended up with a non-square $M$, and that can't happen if you adhere to the sampling theorem.


C/A code (the main civilian channel as of now) has a main lobe of around 2MHz and depending on how much processing you want to do, it is normally sampled in the low MHz range, which you are in. (you could pull it off with 2-5MHz).

Do know that the signal must be mixed down to a lower frequency before you sample it. The original signal is at 1.57GHz.


How to calculate DTFT of cosine function divided by n

The given impulse response is a scaled version of the impulse response of a highpass Hilbert transformer with delay $\tau=N/2$ ($N$ odd) with frequency response
$$H(e^{j\omega})=\begin{cases}-j\,\textrm{sgn}(\omega)\,e^{-j\omega\tau},&\omega_c<\omega<\pi\\j\,\textrm{sgn}(\omega)\,e^{-j\omega\tau},&-\pi<\omega<-\omega_c\\0,&-\omega_c<\omega<\omega_c\end{cases}\tag{1}$$
This can be shown as follows:
$$\begin{align}h[n]&=\frac{1}{2\pi}\int_{-\pi}^{\pi}H(e^{j\omega})e^{jn\omega }d\omega\\&=\frac{1}{2\pi}\left[-j\int_{\omega_c}^{\pi}e^{j\omega(n-\tau)}d\omega+j\int_{-\pi}^{-\omega_c}e^{j\omega(n-\tau)}d\omega\right]\\&=\frac{1}{2\pi (n-\tau)}\left[-e^{j\pi(n-\tau)}+e^{j\omega_c(n-\tau)}+e^{-j\omega_c(n-\tau)}-e^{-j\pi(n-\tau)}\right]\\&=\frac{\cos[\omega_c(n-\tau)]-\cos[\pi(n-\tau)]}{\pi(n-\tau)}\tag{2}\end{align}$$
Since for $\tau=(2k+1)/2$, $k\in\mathbb{Z}$, we have $\cos[\pi(n-\tau)]=0$ we arrive at the following inverse DTFT of $(1)$:
$$h[n]=\frac{\cos[\omega_c(n-\tau)]}{\pi(n-\tau)}\tag{3}$$


How many 8x8 JPEGs exist with distinct DCT given a quality level?

Well, $8\times 8$ is simply 64 individual values.

Now, if you had a quantization matrix $Q$ (dimensions $8\times 8$) with all $1$s, you'd have 8 bit to spend on every entry in your DFT, so that's

$$(256)^{8 \times 8}=({2^8})^{(8 \times 8)} = 2^{512}$$

possible different DFTs.

Now, typically, that's not happening, since quantization matrices restrict the possible values that can end up being non-zero.

The number of possible values is relative easy to figure out: You take the maximum value for each entry, 255, and divide it by the quantization matrix value, round that down; that leads you to the number of possible combinations being

$$ N = \prod_{i = 1, j = 1}^{8,8} \left\lfloor\frac{255}{Q_{i,j}}\right\rfloor$$


Differences in moving window rms due to sampling rate

Shouldn't all the samples times in the set corresponding to the 2 kS
  sampling rate actually be a subset of the sampling times seen in the 5
  kS ?

Not necessarily. It depends on how you do the down sampling. If you use an anti-aliasing lowpass filter, it typically induces time delay and/of phase distortion. If you don't, you probably end up with some residual aliasing, since there is no such thing as a true "band limited" signal.

Could also be rounding errors when you determine the sample locations.


Generally the amount of ringing that you get is a function of the steepness of the filter in the frequency domain, regardless of filter type.
At the same steepness an Elliptic will require a lower order but will have pretty much the same ringing as a Butterworth.

Choosing between a linear phase or minimum phase will change the character but not the extent of the ringing.

No matter how you slice it, a square wave contains a lot of high frequencies and once you filter those out, it's not a square wave any more.


change the frequency component of a spectrum in an audio file

Judging by your question, it sounds like you are developing hearing aid solutions. Given that you are calculating an FFT and are now in the frequency domain, the easiest route forward is to use phase vocoder techniques to shift spectral peaks. Here is the original paper on this technique > http://505606.pbworks.com/f/LaroD99-pvoc.pdf


Relation between filter order and window length for FIR filter

Yes. The number of taps or coefficients in a FIR filter (or window length as you call it) tell the order of the filter. A more general way to see it, to avoid confusions, is to think of the order as the largest delay in the filter. For example, consider the following naive reverb:
$$y[n] = a x[n] + (1 - a) x[n-D]$$
Even though there are only two (visible) coefficients, the order of the filter is $D$ because it’s the largest delay. The window length is $D+1$ but all the middle coefficients are $0$, hence not visible.


Difference between MATLAB hough and my implementation

I don't have the rep to comment on Dima's answer, but the code in the question does convert the $\theta$ to radians already.

This is my implementation of the routine for detecting lines:

[ theta,rho,height,width,Accumulator] = SHTanalyse( m,delta_theta,delta_rho )
%% Get dimensions of input image
[height, width] = size(m);

% Create array for theta
theta = -90:delta_theta:90;
costheta = cosd(theta); % Create lookup tables to optimize computation of cos theta and sin theta
sintheta = sind(theta);

% Find bounds on rho
rhomax = round(sqrt(height^2 + width^2));

% Create array for rho
rho = -rhomax:delta_rho:rhomax;

% Get edgel coordinates
[yEdges,xEdges] =  find(m);

% Initialize the accumulator
Accumulator = zeros(numel(rho),numel(theta));

% Voting process
for k = 1:numel(yEdges) % Loop through edge pixels
    for t = 1:numel(theta)
        % Compute the corresponding rho
        rhvalue = round(xEdges(k).*costheta(t) + yEdges(k).*sintheta(t));
        % Increment accumulator
        Accumulator(rhvalue + round(0.5*numel(rho)) + 1,theta(t) + 90 + 1) = Accumulator(rhvalue + round(0.5*numel(rho)) + 1,theta(t) + 90 + 1) + 1;
    end
end

You can compare its performance with the MATLAB's own Hough function to figure out the discrepancy.


best implementation of a real-time, fixed-point iir filter with constant coefficients

what is fraction saving? can you write a code.so that i can understand more clearly?

Let's call the quantizer operator $\operatorname{Quant}\{\cdot\}$ .  So the output of the quantizer, with $v[n]$ going in, is

$$y[n]  = \operatorname{Quant}\{ v[n] \}$$

which we shall model as an additive error source:

$$y[n]  =  v[n] + q[n]$$

No matter how the quantizer works, we can always derive $q[n]$ in the run-time system with

$$q[n]  \triangleq  y[n] - v[n]  =  \operatorname{Quant}\{ v[n] \}  -  v[n]$$

Even though $q[n]$ is, strictly speaking, a deterministic function of $v[n]$, often, if the signal swings are much larger than the step size, $\Delta$, of the quantizer, we model $q[n]$ as some kinda random number and a "noise source".  It's when the signal gets very small, on par with $\Delta$, that this assumption is no good, and that's when you need to consider adding dither before the quantization operation.  I doubt you will need dither for your situation, but "ya never know".

The simplest quantizer to implement in a fixed-point 2's complement system is simply dropping the less significant bits that are to the right of the rounding point or "binary point":

$$ \operatorname{Quant}\{ v[n] \}  \triangleq  \Delta \left\lfloor \frac{v[n]}{\Delta} \right\rfloor $$

$\lfloor \cdot \rfloor$ is the floor function, the same as the what the C function floor(): returning the most-positive integer that does not exceed the argument to the function.  $\lfloor u \rfloor$ is the sole integer that satisfies:

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

or

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

If $u \ge 0$, often we say that $\lfloor u \rfloor$ is the "integer part" of the real value $u$.  The "fractional part" of $u$ is $u - \lfloor u \rfloor$.  It is clear that the fractional part of $u$ is

$$ 0  \le  u - \lfloor u \rfloor  <  1 $$

So the quantization error is the fractional part remaining after quantization is

$$ q[n]  =  \Delta \left\lfloor \frac{v[n]}{\Delta} \right\rfloor -  v[n] $$

That is the additive noise signal in terms of the input to the quantizer.

$$ -1  <  \frac{q[n]}{\Delta}  =  \left\lfloor \frac{v[n]}{\Delta} \right\rfloor - \frac{v[n]}{\Delta}  \le  0 $$

So we know that the additive quantization error, $q[n]$, has a range of

$$ -\Delta < q[n] \le 0 $$

Now, we assume (that's the only "rule-of-thumb") that essentially the value of $q[n]$ is equally likely in that range so then $q[n]$ is viewed as a random number with uniform p.d.f.:

$$ p_q(\alpha) = \frac{1}{\Delta} \operatorname{rect}\left(\frac{\alpha}{\Delta} + \frac{1}{2}\right) $$

where $\operatorname{rect}(\cdot)$ is the standard unit rectangular function:

$$\operatorname{rect}(u) \triangleq \begin{cases}
1, & \text{if } \ -\frac{1}{2} < u \le \frac{1}{2} \\
0, & \text{otherwise}
\end{cases}$$

The expectation value or mean of $q[n]$ is $-\frac{\Delta}{2}$ (a half-step bias toward the negative direction) and the variance is $\frac{\Delta^2}{12}$.  This is well known about uniform quantization.

In the frequency domain:

$$Y(z)  =  V(z)  +  Q(z)$$

(for frequency response, we substitute $z = e^{j 2 \pi f T}$  where $T$ is the sampling period and in reciprocal units as $f$. Normalized frequency (using the DTFT or "discrete-time Fourier transform") is $\omega \triangleq 2 \pi f T$ and ranges $-\pi < \omega < \pi$.  In the DTFT, if the quantization error really is random, $Q(e^{j \omega})$ would have a Dirac impulse at $\omega=0$ of magnitude $-\frac{\Delta}{2}$ (to represent the DC bias of $-\frac{1}{2}$ LSB because we're always rounding down) and otherwise a flat spectrum from $-\pi$ to $\pi$, that would, because of Parseval's theorem, have area of $\frac{\Delta^2}{12}$.  Maybe there's a $2\pi$ factor in that (which is why i don't usually like angular frequency over "ordinary frequency").

we know that:

$$ \int\limits_{-\infty}^{+\infty} p_q(\alpha) \ d\alpha \ = \ 1 $$

The mean of quantization error (DC component):

$$\begin{align}
\mu_q & \triangleq \overline{q[n]} \\
& \triangleq \lim_{N \to +\infty}\frac{1}{2N+1} \sum\limits_{n=-N}^{N} q[n] \\
& = \int\limits_{-\infty}^{+\infty}  p_q(\alpha) \ \alpha \ d\alpha \\
& = Q(e^{j 0}) \\
& = -\frac{\Delta}{2}
\end{align}$$

The mean square of quantization error:

$$\begin{align}
\overline{\left| q[n] \right|^2} & \triangleq \lim_{N \to +\infty}\frac{1}{2N+1} \sum\limits_{n=-N}^{N} |q[n]|^2 \\
& = \int\limits_{-\infty}^{+\infty}  p_q(\alpha) \ \alpha^2 \ d\alpha \\
& = \frac{1}{2\pi}\int\limits_{-\pi}^{+\pi} \left| Q(e^{j \omega}) \right|^2 d\omega \\
& = \frac{\Delta^2}{3}
\end{align}$$

The variance of quantization error (AC power):

$$\begin{align}
\sigma_q^2 & \triangleq \overline{\left| q[n] - \overline{q[n]} \right|^2} \\
& \triangleq \lim_{N \to +\infty}\frac{1}{2N+1} \sum\limits_{n=-N}^{N} |q[n] - \overline{q[n]}|^2 \\
& = \int\limits_{-\infty}^{+\infty}  p_q(\alpha) \ (\alpha-\mu_q)^2 \ d\alpha \\
& = \overline{\left| q[n] \right|^2} - \left|\overline{q[n]}\right|^2 \\
& = \frac{\Delta^2}{12}
\end{align}$$

Now (to get back to "fraction saving") suppose $v[n]$ is derived from your input, $x[n]$, as so:

$$v[n] = x[n] - q[n-1]$$

So whatever error we added to $v[n]$ by quantizing it, we will subtract, before quantization, from the quantized sample immediately following.

(n.b.: Remember that $q[n] \le 0$ and when we "add" $q[n]$ in the quantization process, we are effectively subtracting whatever bits are to the right of the binary point by zeroing those bits, which is always rounding down.  So when we "subtract" $q[n-1]$ from the quantizer input, $x[n]$, before quantization, we are adding those bits we dropped, extended on the left with zeroed bits, not sign extended.)

$$\begin{align}
y[n] & = \operatorname{Quant}\{ v[n] \} \\
     & = \operatorname{Quant}\{ x[n] - q[n-1] \} \\
     & = \quad x[n] - q[n-1] + q[n] \\
     & = \quad x[n] + (q[n] - q[n-1])  \\
\end{align}$$

or, in the frequency domain:

$$\begin{align}
Y(z) & = X(z) \quad + \quad Q(z) - Q(z) z^{-1} \\
     & = X(z) \quad + \quad (1 - z^{-1}) Q(z) \\
     & =  H_x(z) X(z) \quad + \quad H_q(z) Q(z) \\
\end{align}$$

So the transfer function from input $x[n]$ to output $y[n]$ is still just a wire:

$$ H_x(z)  =  1 $$

But the transfer function from the additive quantization noise $q[n]$ to the output $y[n]$ is

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

So, like FIR filters, there is a pole at the origin $z=0$ (big hairy deel), but there is a zero right at the DC point, $z=1$.  This kills any contribution of quantization error to the output at the frequency of 0.  Even that $-\frac{\Delta}{2}$ bias, from always rounding down, is killed.

The magnitude of the noise to output frequency response is

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

The transfer function, $H_q(\cdot)$, is 0 (or $-\infty$ dB) and when you're up at Nyquist (or $\omega=\pi$) the gain is 2 (or 6 dB).  overall, there is no decrease in quantization noise energy, in fact, because energy is magnitude-squared, there is an increase in total quantization energy, but at frequencies below $\frac{\pi}{2}$, the quantization noise is reduced.  For us audio guys, we like that because none of us hears too well above 11 kHz anyway.

But one thing it helps is that it destroys that limit cycle of IIR filters that gets them stuck on a non-zero DC value, even when the input truly goes to a dead zero.

Code for a biquad (Direct Form I, not Direct Form II which you don't wanna use, especially for fixed-point) that might look like that might be as below. To be clear in sign convention, the transfer function from input to output is:

$$ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}} {1 - a_1 z^{-1} - a_2 z^{-2}} $$

This corresponds to a DF1 difference equation as:

$$ y[n] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] + a_1 y[n-1] + a_2 y[n-2] $$

With quantization it's:

$$ y[n] = \operatorname{Quant}\{ b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] + a_1 y[n-1] + a_2 y[n-2] \} $$

I am sorta assuming that short is 16-bit and long is 32-bit. because, in terms of actual mathematical values, any stable IIR will have  $-2 < a_1 < +2$ and $-1 < a_2 < +1$ ($b_0$, $b_1$, and $b_2$ are relative to each other so we worry less about their fixed range), so let's say that all our coefficients are all scaled up by a factor of $2^{14} = 16384$ (so that the equivalent integer value of the fixed-point number ranges from $-2\cdot2^{14} = -32768 \ <$ a1 $= 2^{14} a_1 < \ +2\cdot2^{14} = +32768$).  so the quantization or rounding point is at bit 14 of the 32 bit number.  The coefficients are also stored as 16-bit short, but we are casting them to long upon loading so that the cast is not continuously necessary in the inner loop.  We want coefficient times state to be a long, even if both multiplier and multiplicand going in are short.

Version 1:

long accumulator;

long a1, a2, b0, b1, b2;   /* 2nd-order IIR coefficients */

short state_x1, state_x2, state_y1, state_y2, state_error;

short x, y;

short get_input(void);

void put_output(short y);

short i, num_samples;
         /* num_samples is the number of samples per block */

/*  load up your coefficients, cast all to long */

/*  load up your states from whereever */

accumulator = (long)state_error;

for (i=0; i<num_samples; i++)
    {
    x = get_input();

    accumulator += b0*x;
    accumulator += b1*state_x1;
    accumulator += b2*state_x2;
    accumulator += a1*state_y1;
    accumulator += a2*state_y2;

    if (accumulator > 0x1FFFFFFF)
       {
       accumulator = 0x1FFFFFFF;     /* clip value */
       }
    if (accumulator < -0x20000000)
       {
       accumulator = -0x20000000;    /* clip value */
       }

    y =  (short)(accumulator>>14);   /* always rounding down */

    state_x2 = state_x1;             /* bump the states over */
    state_x1 = x;
    state_y2 = state_y1;
    state_y1 = y;

    accumulator = accumulator & 0x00003FFF;
    /* keep the fractional bits that you dropped for the next sample
       otherwise clear the accumulator */

    put_output(y);
    }

state_error = (short)accumulator;

/*  save your states back to whereever  */

Now another trick, to speed things up, is to adjust the a1 coefficient by subtracting the scaled fixed-point version of 1 (which, in this case, is $2^{14}$).  Then the code might look like:

Version 2:

long accumulator;

long a1, a2, b0, b1, b2;   /* 2nd-order IIR coefficients */

short state_x1, state_x2, state_y1, state_y2, state_error;

short x, y;

short get_input(void);

void put_output(short y);

short i, num_samples;
         /* num_samples is the number of samples per block */

/*  load up your coefficients, cast all to long */

/*  load up your states from whereever */

a1 = a1 - 0x00004000;    /* adjust a1 coefficient */
/* 0x00004000 is "1" in our scaled-up fixed-point coefficients */

accumulator = (long)state_error + (((long)state_y1)<<14);

for (i=0; i<num_samples; i++)
    {
    x = get_input();

    accumulator += b0*x;
    accumulator += b1*state_x1;
    accumulator += b2*state_x2;
    accumulator += a1*state_y1;
    accumulator += a2*state_y2;

    if (accumulator > 0x1FFFFFFF)
       {
       accumulator = 0x1FFFFFFF;     /* clip value */
       }
    if (accumulator < -0x20000000)
       {
       accumulator = -0x20000000;    /* clip value */
       }

    y = (short)(accumulator>>14);    /* always rounding down */

    state_x2 = state_x1;             /* bump the states over */
    state_x1 = x;
    state_y2 = state_y1;
    state_y1 = y;

    put_output(y);
    }

state_error = (short)(accumulator & 0x00003FFF);

/*  save your states back to whereever  */

Note the missing bit-masking instruction at the bottom of the loop.  So you have the left bits having an extra value of $y[n-1]$ in there.  That's why we have to fix a1 by subtracting 1 (again scaled by $2^{14}$) from it.  That does affect the range of a1 by bumping it over.  Now it's $-2 < a_1 - 1 < +2$ or $-1 < a_1 < +3$ but the range of all of the other coefficients are unchanged.


how to find the frequency of the downsampled signal?

If the term downsampling is used in its correct meaning as anti-aliasing lowpass filtering followed by decimation, then your output should be zero, as the lowpass filter with a cutoff frequency of $\omega = \pi/2$ should remove the input $x[n] = \cos( \frac{2\pi}{3} n) $ ; i.e. , frequency of the cosine is out of the passband of the lowpass filter, and will be blocked by the filter.
However, if you simply mean decimation (or compression) by the term downsampling, then, the output will be: $y[n] = x[2n] = \cos(\frac{2\pi}{3} 2n) = \cos(\frac{4\pi}{3} n) = \cos(\frac{2\pi}{3} n) = x[n]$. This equivalence is due to the aliasing that occurs for omitting the anti-aliasing lowpass filter.
Hence, in this case the frequency of the output is the same as frequency of the input.


An adaptive FIR filter is a FIR filter, that uses some kind of an adaptive algorithm to change the filter weights and reach a desired state. In case of using an LMS algorithm the general update equation is the following:
$$
\mathbf{w}(n+1) = \mathbf{w}(n) + \mu\cdot e(n) \cdot \mathbf{x}(n),
$$
where $\mathbf{w}(n)$, $e(n)$ and $\mathbf{x}(n)$ are respectively filter weights, an error and tap values at $n^{th}$ sample and $\mu$ is an adaptation step size.
That means, that the result of adaptation of a $N$-tap filter will be $N$ new values of its weights.
In your case, you can make a FIR filter $G(z)$ with $N=15$ try to adapt, such that its impulse response will be close to the impulse response of the IIR filter $G_0(z)$. And it should be able to do so starting from zero weights. However, after the convergence it will certainly still have an error.
Useful references on the LMS algorithms are Chapter 9 "Adaptive Filtering" in Statistical Digital Signal Processing and Modeling by M.H. Hayes and a Least_mean_squares_filter entry on Wikipedia.


Is possible to get overshoot on bessel filter?

It's not exactly clear what you mean by overshoot. Let's define it as "when filtering a square wave, the max of of the the output should not exceed the steady state of max of the square wave".

In this case your statement "as i know bessel filter shouldnt have overshoot" is wrong. A bessel filter will indeed overshoot. A few first order filters will not overshoot (butterworth, chebychev, ...) but any filter of order > 1 will overshoot.


While this is by your admission homework (and fairly basic), I'll bite. Recall the definition of the DTFT:

$$
X(\omega) = \sum_{n=0}^{\infty}x[n] e^{-j\omega n}
$$

And recall the definition of the frequency response $H(\omega)$:

$$
H(\omega) = \frac{Y(\omega)}{X(\omega)}
$$

where $x[n]$ is the input to the system and $y[n]$ is its output. Combine these two equations:

$$
\begin{eqnarray*}
H(\omega)X(\omega) &=& Y(\omega) \\
\frac{1 - a^4 e^{-j 4 \omega}}{1 - a^4 e^{-j \omega}}X(\omega) &=& Y(\omega) \\
(1 - a^4 e^{-j 4 \omega})X(\omega) &=& (1 - a^4 e^{-j \omega})Y(\omega) \\
X(\omega) - a^4 e^{-j 4 \omega} X(\omega) &=& Y(\omega) - a^4 e^{-j \omega} Y(\omega)
\end{eqnarray*}
$$

Now, perform the inverse DTFT on both sides of the equation. By definition, $X(\omega)$ and $x[n]$ are a transform pair; likewise for $Y(\omega)$ and $y[n]$. For the other two terms, recall the time-shifting property of the DTFT:

$$
x[n-k] \leftrightarrow e^{-jk\omega} X(\omega)
$$

which can be shown easily from the definition of the DFT. Using this property, the equation inverse transforms to the difference equation specification for the system:

$$
x[n] - a^4 x[n-4] = y[n] - a^4 y[n-1]
$$

$$
y[n] = x[n] - a^4 x[n-4] + a^4 y[n-1]
$$

This is the definition of a recursive filter, which are usually IIR; that is the case for this one. Finding the impulse response is easy; let $x[n] = \delta[n]$ and find that the system output is:

$$
y[n] = a^{4n}u[n] - a^{4(n-4)+4}u[n-4]
$$

$$
y[n] = a^{4n}u[n] - a^{4n-12}u[n-4]
$$

The above is plotted for $a=0.99$. It should be noted that the system is only stable for $|a|\le1$.


Developing a classifier for a periodic music

Very first thing to do: there is literature describing systems doing exactly what you want to do. Look for it!

Your suggested approach - representing the sequence as a feature vector of dimension $n$ where $n$ is the number of onsets or beats - has many flaws:

What if you want to classify sequences which are longer/shorter than your training samples?
What if the onset detector skips an event or detects a spurious event, and the entire vector is shifted left or right?
What about timing information, which is essential in the perception of rhythm!
A single real number computed from the spectrogram is unlikely to represent well the musicological dimension that matters to your problem (which "stroke" or "bol" has been played on the tabla).

Strings/sequence matching methods are more likely to be what you need. What matters to your classification problem is "what follows what" rather than "what happens at the $n$-th position" - because the concept of "$n$-th" position is not very robust, and also very relative.

Here's how I would approach the problem:

Detect onsets on your audio signal, and for each onset, compute a vector of timbral features (MFCCs will be fine). At this point you have a sequence of onset times and a sequence of feature vectors.
On your entire training database, train a vector-quantizer with a small-ish codebook (say 16 or 32 vectors) on the feature vectors. There's a good chance that it will map well to the tabla bols ('dha', 'ti', 're', 'tin'...) - each codeword will correspond to a distinct tabla stroke.
Quantize your data using this vector quantizer. Your sequences will now be represented as a list of onset times and codeword id [(0.0, codeword_10), (0.51, codeword_7), (1.02, codeword_7), (1.27, codeword_16), (1.52, codeword_17), (2.03, codeword_10)...] etc.
Use a classic tactus/tatum detection algorithm (they use inter-onset times as their input) to recover the temporal structure of the rhythm. This will allow you to group events into matras. The example above will become: codeword_10, codeword_7, codeword_7, codeword_16|codeword_17, codeword_10 (note the grouping - the fourth and fifth events are played in sequence within the same matra).
Use any string/sequence modelling technique on the representation above:

Transforming the sequence obtained at step 4 into a sparse vector of n-gram counts is the most flexible method, since it allows you to use any text machine learning technique. The downside is that it needs lots of training data and works well on rather long sequences.
You can also use HMMs. Collect a small corpus of sequences in various taal, train a HMM on them for each taal; use the HMM to compute the likelihood of each of the taal models given the observed sequence.
Another option is to use an edit distance on the sequence obtained at step 4, in conjunction with nearest-neighbour classification. Compute the edit distance between the sequence to classify, and sequences representative of each of the different taals. Pick the taal with the smallest edit distance.

Note that if you have time to invest in manual annotation, steps 2 and 3 can be replaced by a supervised classifier: annotate a database of table sounds with the name or 'bol' of each stroke ('dha', 'dhin', 'ti', 're'...). Train a classifier from audio features (say MFCC) into bols. Use this classifier at step 3. You will now be able to work with a very meaningful representation of your sequences.

A final note: There's a good chance that step 5 can be done with a handful of rules, rather than through machine learning. I know the entire purpose of your project is to do tala recognition with machine-learning; but what is the point if a bit of musicological knowledge can do the same job? What machine learning can do is transcribe the audio signal into a symbolic representation (a sequence of bols) - but further processing of this representation to extract meaningful musicological properties/categories does not necessarily require machine learning techniques...


I've found an estimator that can work for a lot of similar situations, the idea is that you can rewrite the division:
$$ x(t) = z(t) \cdot y(t) \tag1$$
The problem now is to estimate $ z(t) $ that solves the equation.
One can start to find a constant coefficient that minimize the 2-norm. This coefficient will be the linear regression slope without intercept term ($\beta$), hence you can write an approximation of $ z(t) $ that minimize the 2-norm:
$$ x(t) \approx \beta \cdot y(t)\tag2 $$
This estimator is constant, hence it is band limited, and you can estimate the same $ \beta $ from the sampled data, in fact be valid that
$$ x[n] \approx \beta \cdot y[n]\tag3 $$
Obviously this estimation is quite useless but you can extend this approach in a kind of FIR way , first.
You could make a rectangular window and estimate the linear regression slope on $N$ elements.
$$ x[n] \approx \beta(x[(n-N),..,N], y[(n-N),..,n]) \cdot y[n] \tag4$$
I don't like a lot this estimator but should limit aliasing if $ N $ is  large enough.
I instead prefer an IIR like approach, using an exponential window. You can do this using a 0-order RLS algorithm. You can choose a bandwidth choosing the forgetting factor $ \lambda = e^{-\frac{1}{W_s}}$,
where $ W_s $ is the width in samples of the exponential window, that can be expressed in function of cut off pulsation ($\omega_c$) using the relation
$$ W_s = \frac{f_s}{\omega_c} \tag 5$$
where $f_s$ is the sampling frequency.
With RLS you can estimate the slope $\beta$ in a time-variant way, with a limited band in discrete domain, where the bandwidth is limited by $ \omega_c $. Note that in this case (0-order) RLS becomes a scalar problem so it is very easy to implement.


How to perform dechirping of SAR image in MATLAB?

Dechirping is a bit different than matched filtering, but perhaps that not important here.

The code you have looks like it should work if the variables are as you describe. Here is a hypothetical example that range compresses an LFM pulse. The reflected signal from an actual target in the scene is a just a delayed version of the waveform, so it will just be a shifted version of what you get when you range compress the transmitted pulse itself:

Tp = 10e-6;              % Pulse length (s)
Brf = 200e6;             % Bandwidth (Hz)
fs = 1.2 * Brf;          % Sampling frequency (Hz)
Ns = ceil( Tp * fs );    % Number of samples per pulse

% Vector of time samples (s)
t = 1/fs * ( -Ns/2 : Ns/2 - 1 ).';

% The LFM waveform.
lfm = exp( 1j * pi * Brf / Tp * t.^2 );

% Matched filter the pulse.
rngcomp = conv2( flipud( conj( lfm ) ), 1, lfm );

% Display the magnitude of the IPR.
plot( 20*log10( abs( rngcomp ) ) );


Expression for discrete fourier transform of linear ramp

There is indeed a better way; you can derive the analytical expression for the DFT of a ramp. Let's start with the discrete-time Fourier transform (DTFT) of a finite length sequence (i.e., the sequence is zero outside the interval $[0,N-1]$):

$$X(e^{j\omega})=\sum_{n=0}^{N-1}x[n]e^{-jn\omega}\tag{1}$$

We also need the DTFT correspondence

$$nx[n]\Longleftrightarrow j\frac{dX(e^{j\omega})}{d\omega}\tag{2}$$

The DTFT of the constant signal $x[n]=a$, $n\in [0,N-1]$ (and zero otherwise) is

$$X(e^{j\omega})=a\sum_{n=0}^{N-1}e^{-jn\omega}=\begin{cases}aN,&\quad\omega=0\\\displaystyle a\frac{1-e^{-jN\omega}}{1-e^{-j\omega}},&\quad\text{otherwise}\end{cases}\tag{3}$$

From $(2)$ and $(3)$ we get

$$\begin{align}a\sum_{n=0}^{N-1}ne^{-jn\omega}&=j\frac{dX(e^{j\omega})}{d\omega}\\&=-a\frac{Ne^{-jN\omega} (1-e^{-j\omega})-(1-e^{-jN\omega})e^{-j\omega}}{(1-e^{-j\omega})^2},\quad\omega\neq 0\end{align}\tag{4}$$

which is the DTFT of the given sequence. For $\omega=0$ we simply have

$$a\sum_{n=0}^{N-1}n=a\frac{N(N-1)}{2}\tag{5}$$

Now we only need to remember that for finite length sequences, the DFT is a sampled version of the DTFT with $\omega_k=2\pi k/N$. From $(4)$ and $(5)$ we obtain

$$\text{DFT}_N\{an\}[k]=\begin{cases}\displaystyle a\frac{N(N-1)}{2},&\quad k=0\\-a\displaystyle\frac{N}{1-e^{-j2\pi k/N}},&\quad k\in[1,N-1]\end{cases}\tag{6}$$


Parallel processing of large images or large quantities of images is no different than parallel processing for other data types. You will come across the same strengths and limitations.

First of all, there are two major paradigms for parallel problems: Those that are obviously parallelisable and those that are not (or inherently serial).

In an obvious parallel algorithm, the result of a computation can be achieved by the repeated application of one or more functions (the function here as an element of computation) to a dataset. But, in a non-obvious parallel algorithm, the result of the computation cannot be "untangled" and has to be carried out as a single function. There are of course intermediate cases with a task being composed of some parts that can be parallelised in a straightforward way and some tasks carried out serially (unavoidably).

If you are looking for an object in an image through cross correlation, you can apply an image pyramid on it and then feed the different scales to different algorithms that slide the template over the image and try to find interesting points.

If you are trying to derive the "map" of a metric, again, this can be obviously parallelised. For instance, generate an image where each pixel is the standard deviation of the 8 surrounding pixels in the original image.

If you are trying to extract SIFT descriptors, again, this can be parallelised, to a certain extent.

If you are trying to stitch images together, what you cannot parallelise easily is the optimisation process of trying to match the locations of the keypoints. This is because, you are trying to minimise an error that is distributed over many different combinations of keypoints. Therefore, inevitably, even if you tried to do it path-by-patch you would still have to tune the error between patches to make sure that the algorithm has converged to a single optimal solution.

Given what is mentioned in the comments of the question, it might be possible to still converge to a solution as far as the keypoints are concerned. But, re-projecting already large images and then putting them together into one ultra huge image, might still be challenging, not because of lack of parallelism but because of lack of memory.

For more information, please see this link and this link

Hope this helps.

EDIT:

  In the simplest form, I'm attempting to register two gigapixel images A and B

Gigapixel images will require a huge amount of memory if processed in a "conventional" way. This is more or less a fact. The question now is what do we do about it?

I say "more or less" because there are techniques that rely both on software and hardware that could allow in memory processing of such images in the "conventional" way. These are in-memory computing and Shared (or Distributed) Memory. On the former, I forget the exact term right now, but there is a technique which links the actual physical memory of a number of machines together through network links at a very low level. So, suddenly, an algorithm has huge amounts of memory available as if it was "local".

But short of these little beasts the other thing that one can do is process the images in an "unconventional" way. And in this case it really depends on what you want to achieve.

The trick is similar to Sparse Matrices.

Sparse matrices present an interface to the user that is no different than a conventional matrix. You still say I[1200,1500] but internally, the getitem operation of I's data type, that is implied by the use of indexing, does a bit of "magic" with the indices and returns a value. If a matrix is 100x100 and has at most 45 non zero double precision floating point values in it, there is no reason to waste the rest of the space. So, internally, when I[20,32] = 3.1415928 you might store something like 20, 32, [data value] and then later, when 2 * I[20,32], you look up the indices and if they exist in your lookup table you return the value but if not, then you return a 0.0.

In the case of huge images, we are on the other end of the spectrum but still the solution is kind of similar. You create something like a class MappedImage. Instead of loading the image, you have a constructor like MappedImage(pathToDirectory). The contstructor loads the image tiles and does the cross registration but does NOT! actually cross register the images. In other words, all you are after at that point is the transform that takes you from "image space" to "cross registered space".

Once this MappedImage is constructed, you still access it via I[789,145890987] (I here being an object of MappedImage) but the getitem of MappedImage now takes the desired indices, passes them through the transformation matrix, gets the actual locations, "guesses" which image patch contains the pixel values and reads the actual pixel value from that image only.

In this way, you don't have to be keeping all images in memory all the time and you can still use standard memory mapped I/O through the facilities of your language of choice.

The paper provided is primarily focused on doing a high quality registration not necessarily accessing the cross registered image. In a way it does highlight the same problems highlighted earlier and hints at the "conventional" solution which is "work on many small images".

Now, obviously, the next question is "What do you want to do with this 'image'"?

Because, if you want to visualise it, then you can create a MappedImage_Tiler that accesses some I via the interface mentioned earlier and cuts this huge image that you will never be able to see in its entirety into smaller pieces. These can then be sent to an off the shelf tile image viewer and allow you to browse the image like (for example) Google Earth.

But, if you want to work on it and produce some directional derivative or whatever else, then, yes the I[34,2875623897652398756] interface will still work but you will have to lower your expectations on speed. Especially if the goal is to work on "laptop type" hardware as the constraint was in the paper cited.

  I can then map those transformation to higher-resolution levels; however, there is error in this type of mapped registration. I'm hoping to use parallel techniques to register A and B directly at the highest resolution.

This is what I had in mind too: Use the low resolution image derived transform and "project" it to the high resolution images (possibly just through a magnification component (?)).

Hope this helps.


retrieving original data from phase and magnitude of Fourier transform

You need to change how you're calculating the phase from "atan" to "atan2".  The reason is that atan doesn't work in all quadrants.  More on atan2 here.


When a filter is used, we usually design a specific sort of filter: low pass, high pass, or bandpass.  That design assumes an ideal filter that has unity gain in the passband and zero gain in the stop band.

Low order filters to not approximate that "brickwall" very well.  Higher order filters do better.

See this page, which has this diagram showing the ideal lowpass filter and how higher orders approximate it better.


Recording and playback of sound, frequencies question

If your sampling frequency is $f_s=8000$ Hz, your maximum signal frequency is indeed 4000 Hz ($=f_s/2$). If your signal contains frequencies above $=f_s/2$ you would hear the results of aliasing. This means that the original spectrum is folded back into the range $[0,4000]$ Hz. What actually happens is that by sampling with a sampling frequency $f_s$ your spectrum is made periodic with period $f_s$. So if the signal is not band-limited with maximum frequency $f_s/2$, the resulting periodic spectrum will have regions of overlap. E.g., if your signal has frequency components up to 5000 Hz, these frequencies would get folded back to the frequency band between 3000 and 4000 Hz, and consequently lead to distortions in this band. This is the reason why sampling is always preceded by lowpass filtering to attenuate all frequency components that might lead to aliasing (anti-aliasing filter). You can listen to some examples of aliasing here.


Could not construct original matrix using SVD

Don't have enough reputation to comment but you need to use the conjugate transpose in your formula for the result to be correct. So try stftb=U*S*V'; in the last line of code. Note that I removed the . which makes a difference since the matrices you are working with are complex.


Parameter extraction from DFT results, possibly with leakage effects

Interpolation is one method used to estimate parameters of non-periodic-in-FFT-aperture signals.  Parabolic or, better yet, Sinc interpolation will help with estimating both the frequency and the amplitude peak B.  For phase interpolation, I find it easier to use something like a pre or post fftshift to move the phase 0 reference to the middle of your data, then interpolating an estimated phase from the fft, then offsetting the phase back relative to wherever using the frequency estimate.


Reliably matching a degraded image with the original out of dozens of candidates

Your question has many solutions which can be good or bad, so I am giving the solution I came up with, hope it is a good one.
Since your images are only different in resolution and color I think the best way to compare them is to make their size the same and then to get rid of color and shading differences only compare their edges not the images themselves. This way, you can use comparison criteria that you've tried now independent of color and shading variations. I simply used sum of differences between two images in the code below and for the images you provided seems it works.

Through this code I try to give you an idea idea on how it is working:

I1=rgb2gray(imread('Low.png'));
I2=rgb2gray(imread('High1.png'));
I3=rgb2gray(imread('High2.png'));

[h,v]=size(I1);
I2=imresize(I2,[h,v]);
I3=imresize(I3,[h,v]);

d1=edge(I1,'log');
d2=edge(I2,'log');
d3=edge(I3,'log');

score1=sum(sum(abs(d1-d2)))
score2=sum(sum(abs(d1-d3)))

Some low pass filtering before calculating the scores might help.


Goodness of fit for complex valued curves (i.e. frequency responses in frequency domain)

You could use the $p$-norm of the complex error $E(k)=Y_{data}(k)-Y_{syn}(k)$:

$$\epsilon_p=\sqrt[p]{\sum_{k=1}^K|E(k)|^p}\tag{1}$$

The most common values for $p$ are $p=2$ (Euclidean norm, $l_2$-norm) and $p=\infty$ (maximum norm or Chebyshev norm):

$$\epsilon_2=\sqrt{\sum_{k=1}^K|E(k)|^2}\tag{2}$$
$$\epsilon_{\infty}=\max_{1\le k\le K}\{|E(k)|\}\tag{3}$$

It eventually depends on your application which error measure best reflects the quality of the approximation. Note that for the maximum norm $(3)$, a single index $k$ determines the value of the error measure $\epsilon_{\infty}$, whereas the error measure $\epsilon_2$ in $(2)$ depends on all values of the complex error $E(k)$, and it is proportional to the square root of the mean squared error.


Welch's overlapped method for amplitude by frequency estimation

I've found this:
Implementing Welch's method for Power Spectral Density

That seems to answer the question. However, the author encountered a problem, which I believe may be addressed by modifying the code that was identified as problematic in the thread by:

dataA(istart:iend,j) = data1(istart:iend,j).*window;
dataB(istart:iend,j) = data2(istart:iend,j).*window;


Is Hilbert transform not defined for complex signals?

The Hilbert transform can be applied to complex functions of a real variable. E.g., the Hilbert transform of the complex exponential $e^{j\omega_0t}$, $\omega_0>0$, is given by

$$\mathcal{H}\{e^{j\omega_0t}\}=-je^{j\omega_0t},\qquad\omega_0>0$$

The problem you encounter has to do with Matlab's implementation of the function hilbert.m. It is designed for real-valued input sequences and it will ignore any imaginary part. Note that despite its name this function does not simply return the Hilbert transform of the input vector, but it computes the corresponding analytic signal, i.e. it returns a complex vector, the real part of which is equal to the input vector, and the imaginary part of which is the Hilbert transform of the input vector.

So if for whatever reason you want to compute the Hilbert transform of a complex vector, you need to do the following:

x = ...    % some complex vector
xr = real(x);
xi = imag(x);
xr_ = imag(hilbert(xr));    % Hilbert transform of real part
xi_ = imag(hilbert(xi));    % Hilbert transform of imaginary part
x_ = xr_ + 1i * xi_;        % Hilbert transform of complex vector x


I think the reason you are getting confused is that you are only considering the convolution of two signals, whereas convolution is not usually between signals but between a signal and an impulse response.

The best way to visualise it is probably with a filter.  For example, in

we can see that the input comes in, and is multiplied by $b_0$ to give its contribution to the output $y$.  Then in the next time step, that input value is multiplied by $b_1$ to give its contribution to the output, and so forth.  In this case your $k$ variable specifies which time in the output sequence you want.  The output signal is given by $y[n] = (x * b)[n]$.

For example, consider the case of $n=4$ (that's the same as $k=4$ in your example).  That means we're looking at the 4th time step (since time = 0).  So the value at the input is $x[4]$, the fourth input value; the value to the right of the first $z^{-1}$ block is the previous input value $x[3]$; the value to the right of the next $z^{-1}$ block is $x[2]$; and so forth.  The output $y[4]$ will be given by
$$y[4] = x[4]b_0 + x[3]b_1 + x[2]b_2 + x[1]b_3 + \cdots + x[-3]b_7 + \cdots$$
 which you should be able to see is the discrete convolution sum, and includes the product of $x[-3]$ with $b[7]$.

Any linear filter or transfer function can be implemented as a convolution, where the output is given by the convolution of the input with the system impulse response.  Try reading up on convolution in a text book, and look particularly at the sliding-strip method of calculating it, which might help with your visualisation.


Difference between Nyquist rate and Nyquist frequency?

Harry Nyquist invented/discovered/proved a lot of things; it can be hard to keep track of them all. The three most important for signal processing and communications are probably these:

If you sample a (real) signal $s(t)$ at $f_s>2B$ samples per second, then $s(t)$ can be reconstructed from its samples, where $B$ is the bandwidth of $s(t)$. The lower bound $2B$ is often called the Nyquist rate.
If you sample a signal $s(t)$ at $f_s$ samples per second, then $f_N=f_s/2$ is frequently called the Nyquist frequency (also folding frequency). All frequencies in $s(t)$ higher than $f_N$ will be aliased back to a frequency between 0 and $f_N$.
Assume you have a baseband channel of bandwidth $B$ (that is, signals with bandwidth $B$ or less are not distorted by the channel). This channel allows transmission of a maximum of $2B$ orthogonal pulses per second. $2B$ in this context is also called the Nyquist rate.


Filtering out a specific Gaussian from a signal using Fourier concepts

You can try performing a matched filtering with the gaussian of interest. See the below python code for an example. It is not exactly a work in the Foruer domain, (but can be seen as such) and remove thes noise just leaves the Gaussians. You need to adapt the sigma of your gaussians (i just did some rough trial-and-error here)

signal = np.loadtxt("C:/local/mmatthe/work/sheet1.csv")
t = np.arange(-signal.shape[0]/2, signal.shape[0]/2)
t2 = np.arange(-20, 20)
sigma2 = 10
gauss = lambda t: 1/np.sqrt(2*np.pi)*np.exp(-(t)**2/(2*sigma2))
gauss_samples = gauss(t2)
gauss_samples /= np.sqrt(sum(gauss_samples**2))

plt.figure(figsize=(10, 5))
plt.subplot(1,2,1)
plt.plot(t, signal)
plt.plot(t, gauss(t))

plt.subplot(1,2,2)
plt.plot(np.convolve(gauss_samples, signal));


I think you made a mistake between the power spectrum and the Fourier Transform.
This is the right form of the Convolution Theory (Nothing is squared):

$$ \lvert Y(\omega)\rvert = \lvert H(\omega) X(\omega)\rvert $$

Try this and it will work for you on MATLAB.

MATLAB Code

% Calculating the Magnitude of the Filtered Signal

vFilterCoeff = [1; 2; 3; 4; 5];
vInputSignal = randn(1000, 1);

vOutputSginal = conv(vInputSignal, vFilterCoeff);

fftNumPoints = 1024;

vOutputSignalDftMagnitude = abs(fft(vOutputSginal, fftNumPoints));
vFilteredSignalDftMagnitude = abs(fft(vFilterCoeff, fftNumPoints) .* fft(vInputSignal, fftNumPoints));

hFigure = figure();
hAxes   = axes();
plot([vOutputSignalDftMagnitude, vFilteredSignalDftMagnitude]);

The result is:

As can be seen there's a full equality.


Why is my signal attenuated after applying a HighPass-filter?

The correct answer to this question has been provided by: Stanley Pawlukiewicz but has since been removed.

  Your filter attenuates the high frequencies of the total timeseries,
  signal and noise.

  Your results are expected.

  You really can’t perfectly remove just the noise.

  When the noise and signal overlap in time and frequency, you have to
  make a tradeoff between attenuating signal and attenuating noise, by
  adjusting the cutoff frequency.


How does the 'phase' of a cosine wave become 'pi' when it transitions from positive to negative?

The statement in the video is, at best, misleading. Normally, when people speak of the "phase" of a sinusoidal function $\cos \omega$, they are talking about the value of $\omega$. For instance, when we say that a wave undergoes a "phase shift", we mean that the value of $\omega$ has changed. More information on this definition of the word "phase" is at https://en.wikipedia.org/wiki/Phase_(waves).

As noted by Jazzmaniac in comments, we could also try to interpret the statement by regarding the value of $\cos \omega$ as a complex number $z$. In polar representation we can write $z = |z| e^{i \phi}$ where $\phi$ is the phase of the number $z$. This use of the word "phase" is misleading in this context. For instance, applying a "phase shift" would mean that we multiply $\cos \omega$ by some complex number $e^{i \phi'}$. That certainly wouldn't implement a phase shift of the sinusoidal wave, though, in the usual sense.


I guess that the missing area is due to a ring. A solution would be to use a PCA in order to compute the finger orientation. Then you apply a closing (mathematica morphology) with an segment oriented segment (same orientation as the finger).

It might modify a little bit the gap between the fingers where they are connected to the hand.


Why is carrier wave affected by percentage modulation?

A conventional AM signal can be written as

$$x_{AM}(t)=A_c\big(1+\alpha\; m(t)\big)\cos(\omega_ct)\tag{1}$$

where $\alpha$ is the modulation index and $m(t)$ is the message signal. In Eq. $(1)$ it is assumed that the message signal $m(t)$ is normalized, i.e., $\max_t|m(t)|= 1$. The average power of $x_{AM}(t)$ is

$$\overline{x^2_{AM}(t)}=\frac{A_c^2}{2}\left(1+\alpha^2\overline{m^2(t)}\right)\tag{2}$$

If you assume now that the power $(2)$ is fixed to some value $P$, then the total available power must divided between the carrier term and the message-bearing term in $(2)$. Equating $(2)$ with the available power $P$, and expressing the carrier term $A_c^2/2$ in terms of the modulation index and the average power of the message signal gives

$$\frac{A_c^2}{2}=\frac{P}{1+\alpha^2\overline{m^2(t)}}\tag{3}$$

In the figure in your question it looks like they assumed $\overline{m^2(t)}=\frac12$, which for a maximum modulation index $\alpha=1$ results in

$$\frac{A_c^2}{2}=\frac23 P\tag{4}$$

This value corresponds to the right-most point of the curve labeled "carrier".


What does it mean by "Adding Gaussian Noise to a 2D Image"?

You can find Gaussian noise generators in many image processing softwares.At each pixel, you add a realization from such random noise generators. You can look for more details at How to Generate White Gaussian Noise.

On StackExchange, additional sources of information:

Why Gaussian noise is usually used?
Why is Gaussian noise called so?
What's gaussian noise look like?
The model of noise in Image denoising problem


Detecting and isolating part of an image

Since the Extraneous part is connected with leg section via a small thread like thing, Try converting it to binary Image, and then eroding the image to get rid of smaller objects. That should disconnect your both objects. Then choose the object with biggest area.

Helpful functions - im2bw, imerode, rogionprops.


How does taking the absolute value of a complex signal reflect in the frequency domain?

Absolute value is not a linear operation, so no, there is not a consistent dual in the frequency plane.

Having said that, though, it is somewhat similar to calculating the signal power, which involves multiplying $x[n]$ with it's complex conjugate.  Multiplication in the time domain is equivalent to convolution in the frequency domain.  The result in the time domain is real, so the result in the frequency domain is even (i.e. symmetric around 0 Hz).

Absolute value will probably produce a similar effect to calculating the power since it is the square root of the power.  It is also real so it will be even in the frequency plane.


Is there a specific reason you need to use pulstran() ? Since you have an explicit expression for your signal, you can just sample the sine wave at the specified times as follows:

function x = ssin(f0,Fs,N)

% f0 = sine wave frequency [Hz]
% Fs = sampling frequency [Hz]
% N = number of samples

dt = 1/Fs;       % Sampling period [s]
T = N*dt;        % Record window duration [s]
t = 0:dt:(T-dt); % Time vector for sampling [s]

% generate samples at the specified times
x = sin(2*pi*f0*t); % [signal units]

stem(t,x)

end

% run the function
x = ssin(1,20,60)


Proof that spectrum analyzer measures mathematical Power Spectral Density

I am assuming $X(t)$ is a passband signal having finite support (non zero value for only a finite range of frequencies within spectrum analyzers range).

In your equation (1) what you have missed is $f_{LO}$ to $f_{IF}$ conversion. So $W_{\Delta t}$ should be $W_{IF}(t) = e^{j2\pi (f_{LO}-f_{IF})t}W_{\Delta t}(t)$. When you convolve $X(t)$ and $W_{IF}(t)$, in frequency domain you will get the magnitude as $|(X(f - f_{\Delta}) * W_{IF}(f - f_{\Delta})|^2$. Here $f_{\Delta}$ is the frequency shift needed to bring $X$ to IF and $W$ also includes the windowing parameter necessary to remove discontinuities around $\pm\Delta t/2$.

In (2), you are first applying a mixer-filter but then comparing (1) - a frequency domain quantity with (2) a time domain quantity. Let us say, the quantity in (2), you have applied the windowing after mixing to IF. ie, $X_{IF}(t) = (X(t)*F_{IF}(t)) \times W_{\Delta t}(t)$.

If you take fourier transform of above, you will get

$|(X(f - f_{\Delta }) F_{IF}(f - f_{\Delta})) * W_{\Delta t}(f - f_{\Delta })|^2$. But for the frequencies of interest (remember the spectrum analyzer has only finite bandwidth, so it will filter out frequencies outside it) the value of $F_{IF}(f) = 1$. So finally it will be $|X(f - f_{\Delta}) * W_{\Delta t}(f - f_{\Delta})|^2$ which is the same as we got above. Your intuition is right about (2). It is indeed convolution in time domain. But if transform it back to frequency domain, you can see it is same as (1). There are lot of things happening in the background. We have sampling and before sampling we have the Low Pass filter. All this will eliminate the unwanted frequencies outside the band of interest. Hence we can 'forget' the effect of filter $F_{IF}$


Single-pole IIR low-pass filter - which is the correct formula for the decay coefficient?

The given single-pole IIR filter is also called exponentially weighted moving average (EWMA) filter, and it is defined by the following difference equation:

$$y[n]=\alpha x[n]+(1-\alpha)y[n-1],\qquad 0<\alpha<1\tag{1}$$

Its transfer function is

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

The exact formula for the required value of $\alpha$ that results in a desired $3\;\textrm{dB}$ cut-off frequency $\omega_c$ was derived in this answer:

$$\alpha=-y+\sqrt{y^2+2y},\qquad y=1-\cos(\omega_c)\tag{3}$$

Even though it should be easy enough to compute $\alpha$ from $(3)$, there are several approximative formulas floating around the internet. One of them is

$$\alpha\approx 1-e^{-\omega_c}\tag{4}$$

In this answer I've explained how this formula is derived (namely, via the impulse invariant transformation of the corresponding continuous-time low pass filter). This answer compares the approximation $(4)$ with the exact formula, and it is shown that $(4)$ is only useful for relatively small cut-off frequencies (of course, small compared to the sampling frequency).

In the wikipedia link given in the question, there is yet another approximative formula:

$$\alpha\approx\frac{\omega_c}{1+\omega_c}\tag{5}$$

[Note that in all formulas of this answer $\omega_c$ is normalized by the sampling frequency.] This approximation is also derived from discretizing the corresponding analog transfer function, this time not via the impulse invariant method, but by replacing the derivative by a backward difference:

$$\frac{dy(t)}{dt}{\huge |}_{t=nT}\approx\frac{y(nT)-y((n-1)T)}{T}\tag{6}$$

This is equivalent to replacing $s$ by $(1-z^{-1})/T$ in the continuous-time transfer function

$$H(s)=\frac{1}{1+s\tau}\tag{7}$$

which results in

$$H(z)=\frac{1}{1+(1-z^{-1})\tau/T}=\frac{\frac{1}{1+\tau/T}}{1-\frac{\tau/T}{1+\tau/T}z^{-1}}\tag{8}$$

Comparing $(8)$ to $(2)$ we see that

$$\alpha=\frac{1}{1+\tau/T}\tag{9}$$

Since the (continuous-time) $3\;\textrm{dB}$ cut-off frequency is $\Omega_c=1/\tau$, we obtain from $(9)$

$$\alpha=\frac{\Omega_cT}{1+\Omega_cT}\tag{10}$$

Equating the discrete-time cut-off frequency $\omega_c$ with $\Omega_cT$ in $(10)$ gives the approximation $(5)$.

The figure below shows the actually achieved cut-off frequency for a given desired cut-off frequency for the two approximations $(4)$ and $(5)$. Clearly, both approximations become useless for larger cut-off frequencies, and I would suggest that approximation $(5)$ is generally useless.


OBPSK in the way that you describe it is also called $\pi/2$-BPSK meaning that the signal constellation is different from bit to bit, rotating by $\pi/2$ after it is used. There are two variations of this format.

Use $\pm \cos(\omega_0 t)$ to transmit one bit and
$\mp \sin(\omega_0 t)$ to transmit the next bit. Repeat this over and over again.  Note the $\mp$ before the $\sin$ carrier.
This method rotates the constellation by $\pi/2$ counterclockwise in going from $\pm \cos(\omega_0 t)$ to $\mp \sin(\omega_0 t)$, and clockwise when going from $\mp \sin(\omega_0 t)$ to $\pm \cos(\omega_0 t)$.
Use

$\pm \cos(\omega_0 t)$ to transmit the first bit
$\mp \sin(\omega_0 t)$ to transmit the second bit
$\mp \cos(\omega_0 t)$ to transmit the third bit
$\pm \sin(\omega_0 t)$ to transmit the fourth bit,
and repeat this over and over. Note the $\mp$ and $\pm$ signs. This
corresponds to the constellation rotating by $\pi/2$ counterclockwise after each and every bit transmission.

Neither form of OBPSK signal can be demodulated very well
using a OQPSK demodulator, though it would be possible to do
so. In OQPSK,
even-numbered bits are transmitted using $\pm\cos(\omega_0t)$ and
odd-numbered bits using $\mp\sin(\omega_0 t)$ but the important
difference is that the $2n$-th bit is transmitted during the
time interval
$[2nT, (2n+2)T)$ while the $(2n+1)$-th bit is transmitted during
the time interval $[(2n+1)T, (2n+3)T)$, i.e. in overlapping but offset time intervals.
At the demodulator, the matched filters/correlators in the I branch
and Q branch thus integrate/correlate over offset time intervals,
and the demodulated data bits come out of the OQPSK demodulator alternately on the I and Q lines, one every T seconds.
On the other hand, in
OBPSK, the $2n$-th bit is transmitted during $[2nT, (2n+1)T)$ and
the $(2n+1)$-th bit during $[(2n+1)T, (2n+2)T)$. If demodulated with
a plain vanilla OQPSK demodulator, during $[(2n+1)T, (2n+2)T)$ the matched filter/correlator input in the I branch would be just noise
(and similarly for the Q branch).
This will decrease the effective SNR. So, Yes, a standard OQPSK
receiver could be used to demodulate an OBPSK signal, but
there are better ways of demodulating an OBPSK signal than just
processing it through an OQPSK receiver.


Unstable frequency/phase when using external reference for Ettus B200

You're looking only at the real part. Obviously, if the phase of transmitter and receiver aren't locked (like here), that will mean the sent oscillation will be more or less in-phase.


What do I lose when I set overlap rate to a really high value in STFT

If you don’t care about computational costs, you can start a window at each sample (e.g. 100% - 1 sample overlap).  It’s well into diminishing returns, but phase vocoder estimation methods work slightly better with greater overlap.


How to compute autocorrelation of signal defined by difference equations?

There are different approaches to compute the requested auto-correlations but I would like to provide the simplest that suits to your case provided that.

The difference equation is an LCCDE type indicating an LTI system.
The input to the system $v[n]$ is WSS (wide sense stationary) random process.

Both of which I assume to be observed in your question. Then I would use the following well known relation between the input and output auto-correlations of an (real) LTI system driven by (real) WSS input $v[n]$:

$$ r_{dd}[k] = h[k] \star h[-k] \star r_{vv}[k] $$

where the auto-correlation of WSS $v[n]$ is $r_{vv}[k] = \sigma_v^2 \delta[k] = \delta[k]$ for an uncorrelated, zero mean, unit variance process. And $h[n]$ is the impulse-response of the LTI system signified by the LCCDE.

Then the auto-correlation sequence of the WSS output $d[n]$ is
$$ r_{dd}[k] = h[k] \star h[-k] = \sum_{m=-\infty}^{\infty} h[m]h[-(k-m)]~~~,~~~\text{ for } k = 0,\pm 1, \pm 2,... $$

So assuming that your LCCDE :
$$d[n] - \alpha  ~d[n-1] = v[n] ~~~, ~~~\text{ with } |\alpha| < 1 $$

signifies a causal LTI system, then its impulse response is:

$$ h[n] = (-\alpha)^n u[n] $$

and yields the auto-correlations to be:
$$ r_{dd}[k] = \sum_{m=-\infty}^{\infty} (-\alpha)^m u[m] (-\alpha)^{m-k} u[m-k]~~~,~~~\text{ for } k = 0,\pm 1, \pm 2,... $$

I hope you can proceed the rest.


What does "convergence" mean for Thinning?

The thinning operation become idempotent after a number of iterations, which means the result does not change anymore - in other word, the algorithm have converged to a solution.


"General" DSP (and signal processing in general) tends to deal much more with problems that can be solved with linear filters, or at least a preponderance of linear filtering, and there's a lot more content involving 1-D signals -- so you tend to see a lot of emphasis on frequency-domain analysis and FFT (and even some folks who can't seem to distinguish between typing 'fft' in Matlab and understanding DSP).

Image processing (video processing in particular) tends to deal with problems that must be solved with lots of nonlinear processing -- image segmentation, finding edges, etc.  Even super-duper neural networks have significant nonlinearities.  While there are linear operations (i.e., convolution, cross-correlation) these tend to be processor-intensive, and are used sparingly.

One can go through an entire career in DSP and not do image processing, and visa-versa (I'm a case in point, with around 30 years experience in DSP and digital control systems, yet with very little image processing under my belt).

About the only thing that is really common, besides the fact that you're taking in something you could call a "signal" and you're doing something to it you could call "processing", is that "plain old" DSP and image processing both reward the practitioner who's willing to do the hard math (or who's lucky enough to find the hard math easy).

It also helps to be able to make intuitive leaps to potential solutions, and then do the hard math to back-fill to determine that your solution will actually work, to refine it, and to be on a firm mathematical foundation for your next intuitive leap.


Optimal $ n $ -th Order IIR /AR Approximation of a Moving Average Filter

(Update: I just realized that the first part of this covering CIC structures is basically what Hilmar has already answered-- I'll leave this up since it offers more graphics and details in case that his helpful to anyone but it is indeed the same answer)

This may not be optimum and although a feedback structure is involved it is strictly not IIR, but want to point out the efficiency of the cascade integrator-comb structures for implementing boxcar moving average filters- and importantly how the results match exactly and not an approximation. So a 600 tap moving average filter can be implemented with simply one subtraction after a 601 sample delay followed by an accumulator.

These structures are quite popular when resampling is involved (often after a moving average the higher sampling rate of the input is no longer required and a lower output sampling rate would be even more efficient).

The equivalency is easily proven from the Z-transform for each of these structures. Specific to the example structures in the graphic:

The Z-transform for the accumulator is:

$$\frac{1}{1-z^{-1}}$$

The Z-transform for the comb (delay and subtract) filter is:

$$1-z^{-4}$$

And the Z-Transform for the moving average filter is:

$$1+z^{-1}+z^{-2}+z^{-3}$$

Combining confirms the equivalence:

$$1+z^{-1}+z^{-2}+z^{-3} = \frac{1}{1-z^{-1}}(1-z^{-4})$$
$$(1-z^{-1})(1+z^{-1}+z^{-2}+z^{-3})= (1-z^{-4})$$

Another favorite of mine as a true IIR alternative to a moving average filter is the exponentially weighted moving average filter, where assuming sufficient bit precision, very tight (in other words effectively long delay) filtering can be achieved by using an alpha that approaches 1. This is demonstrated in the following graphics:

In fixed precision structures the order of the two operations shown in the block diagram of the exponentially weighted moving average filter is very important! In most cases when dealing with fixed precision it is best to allow digital filters to grow the signal first and do the scaling down after the accumulation. This also explains the motivation for extended precision accumulators- the accumulator shown needs to be designed with sufficient precision to prevent any overflow for the given $\alpha$ value chosen (the Infinite Summation identity $\sum_{n=0}^\infty \alpha^n =\frac{1}{1-\alpha}$ is useful in this case to determine worst case maximum value; resulting in $gain = \frac{1}{1-0.99}$ or 100 in the example shown, so in this case we would need to add $log2(100) = 7$ bits in the accumulator). Also clear from the summation identity is how the scaling term itself was determined such that the gain of the structure is 1.
The order of these operations is no different than the reason to put an amplifier prior to an attenuator in an analog system to get the best noise figure, while paying attention to clipping and maximum input signal, for those familiar with that.


Is any time-domain vector stored on a digital machine always intrinsically causal?

Causality is not so much a characteristic of a signal as it is a characteristic of a system.  For example, a non-causal system can have an output at time $t$ which depends on the input at time $t+1$.  When thinking in terms of time, a non-causal system breaks our intuition because it has to "see the future" in order to operate.

Let's say that I want to create an audio effect that is a sort of "reverse echo."  In other words, I want to hear an echo of the sound before the sound event actually occurs.  This would be an example of non-causal processing because at any given time, the output depends on input which has not yet occurred. This is not a problem in the case of a recorded audio signal because we already have all of the time samples available to us, so we can "look ahead" and use this "future information" right now.

But what if I wanted to implement this "reverse echo" effect in a live performance?  I obviously can't "look ahead" to grab samples out of the future, but I can wait until all of the samples I need are available to me, and then apply the processing.  This would produce the exact same output signal as if I had done the non-causal processing mentioned above, but with one significant difference: my output would be delayed.

Given a signal $y$ which depends on input $x$ as

$$y[n] = a_2 x[n-1] + a_1 x[n] + a_0 x[n+1]$$

we can make a causal version of $y$, which we will call $y^\prime$, by simply delaying $y$ by one time step:

$$y^\prime[n] = y[n-1] = a_2 x[n-2] + a_1 x[n-1] + a_0 x[n]$$

We have not altered the signal in any way, besides to shift its time indexing. Any shift in phase can be seen as a direct result of the delay, and cannot be reversed by multiplying some phase term.


Poisson noise and curve fitting - denoise first?

If the example images you've given are at all representative of your application, you may want to consider thinking about the problem a little differently. Instead of thinking of the image as "corrupted by Poisson noise", think of the observed data as a limited number of photons sampled at random from the latent image intensity map. The photon counts you get are providing incomplete information about the latent image, not corrupt or noisy information.

From this perspective, if you are truly fitting a curve or formula of few parameters to your data, this is a classical parametric statistical estimation problem (lots of samples, few parameters). Such problems have a good estimation theory which guarantees you will accurately learn the fit parameters without needing too many sample points. If you apply a denoising step before parametric fitting, the theoretical guarantees will no longer apply because the error in the denoised image is some complicated composite of residual fit variance and bias from the denoiser's implicit signal model.

So there's a theoretical justification for fitting the function straight to the data. Just make sure that you use the correct Poissonian likelihood function when you do so, or use the EM algorithm if you don't want to worry about such things.

EDIT: actually you don't even need Poissonian likelihood here (unless your data is quantized/binned). Just take the approach described in http://en.wikipedia.org/wiki/Maximum_likelihood. Here the $x_1, x_2, ... x_n$ are a collection of $n$ vectors giving the 2-D coordinates of your photon observations, and $\theta$ is the parameter vector for the function you want to fit. In the case of a Gaussian you have $\theta = (\mu,\Sigma)$ where $\mu$ is the mean (centroid in 2D) and $\Sigma$ is the covariance matrix. Then in the notation of the wiki article, your likelihood is

$$ \prod_{i=1}^n f(x_i | \theta) = C \prod_{i=1}^n \exp( -\tfrac{1}{2} (x_i - \mu)^T \Sigma^{-1} (x_i - \mu)) $$

where $C$ is some constant with square roots of pi in it and such :) If your data is quantized/binned then you consider the count value in each bin as a Poisson random variable and go from there. Weirdly, the Poisson and the $n$ observation approaches are nearly equivalent from an optimization standpoint.

I don't have a great reference for you offhand because I learned this stuff in lab classes or wherever years ago. But the wiki article and refs therein seem good. This paper seems ok too: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.74.671&rep=rep1&type=pdf


Difference between slow, medium, fast AGC

I don't believe there is much more to it than that, especially noting that you are aware not to track so fast that you track out amplitude information for modulated signals where the amplitude is used as part of the modulation (so more so for M-QAM than M-PSK). You can also design AGC's with different rise/fall times so can have "fast-attack" and "slow-decay" for example.

A reason you may want to have a dynamic AGC is distinguishing operation between acquisition and tracking. During acquisition it would be likely to use a much faster AGC at the expense of noise performance for the benefit of acquiring the signal faster. Then once acquired the loops (AGC, carrier tracking, symbol tracking) can run slower at the benefit of better noise performance but fast enough to track all expected dynamics.