Finding set of orthogonal basis functions for composite signals I've been given a list of 5 composite signals, where each is composed of 10 sinusoids of different frequencies. For instance, the first composite signal $S_1$ is given by
$$
S_1 = \sum_{i=1}^{10} A_i \sin (2 \pi f_i t- \phi_i)
$$
where $A_i, f_i$ and $\phi_i$ are known values. The other 4 composite signals are similar to $S_1$. In total, I have 50 different values for $A_i, f_i$ and $\phi_i$ (10 for each $S_i$).
My goal is to find the orthogonal basis functions for $\{S_1, S_2, \cdots, S_5\}$. I've decomposed them into $\sin$ and $\cos$:
$$
\sum_{i=1}^{10} A_i \sin (2 \pi f_i t- \phi_i) = \sum_{i=1}^{10} A_i [\sin(2\pi f_i t)\cos(\phi_i) - \cos(2\pi f_i t)\sin(\phi_i)]
$$
and from here I can see that $A_i, \cos(\phi_i), \sin(\phi_i)$ are constants, and a spanning set would be
$$
\bigcup_{k=1}^{50} \{\sin(2\pi f_k t), \cos(2\pi f_k t) \}
$$
so I end up with a 100-element set. I presume these are my orthogonal basis functions, since each $f_k$ is unique. Can I check if my train of thought is right?

Your $S_1,\ldots, S_5$ already are your orthogonal bases, exactly the way your first equation writes them down.
To check that: write down the definition of your inner product, and insert, for example, $S_a$ and $S_b$. You'll see that $\left<S_a,S_b\right>=0 \iff a\ne b$.
This follows quite directly from the linearity that a (real) inner product needs to have:
\begin{align}
\Big<\sum_i A_{a,i}\sin(2\pi f_{a,i}t -\phi_{a,i}),
\sum_k A_{b,k}\sin(2\pi f_{b,k}t -\phi_{b,k})
\Big>&\\
= \sum_i
\Big<
A_{a,i}\sin(2\pi f_{a,i}t -\phi_{a,i}),
\sum_k A_{b,k}\sin(2\pi f_{b,k}t -\phi_{b,k})
\Big>&\\
=\sum_i\sum_k
\Big<
A_{a,i}\sin(2\pi f_{a,i}t -\phi_{a,i}),
 A_{b,k}\sin(2\pi f_{b,k}t -\phi_{b,k})
\Big>
\end{align}
and sines of different frequencies are orthogonal.


Pulse Compression(Chirp) SNR Gain I've created both single tone signal and a LFM chirp signal with same pulse duration and cross correlated both of them with their own(matched filter). Peak points of both results were same as expected which is the total energy of input pulse. I've read about SNR gain with signal compression in radar equations and chirp is an example of compressed signal. Where does this SNR gain come from? Does it simply mean matched filter? If so pulse compression is no different than a single tone as far as I understand(range resolution is the only difference).

Edit:
Equation 7.4 states that output peak SNR to input SNR is proportional to both bandwidth and pulse length. It looks at the output to input SNR ratio which is affected by the input noise because of bandwidth used at the input. Does this mean a single tone signal has a better input SNR because of less bandwidth? Then what is the paper trying to say by comparing output SNR to input SNR if peak point is the total energy of my pulses on both compressed and single tone after a matched filter?

Edit2: After some thinking I've concluded it as follows: No matter what the signal modulation or bandwidth, peak SNR at matched filter output is same for input pulses of same energy. However, input noise power is proportional to the bandwidth of signal which gets filtered on matched filter and that's what we call matched filter gain. Correct me if I'm wrong.

It's often said that pulse compression gives you a gain proportional to the time-bandwidth product (otherwise known as the pulse compression ratio, or $PCR$). This is a really misleading statement, and it had me confused enough to sit down and think about it for awhile. I thought I'd share some of my findings that I pieced together from both reading the internet and from my own DSP simulations.

The claim

Before looking into this, I had assumed that, if I take two pulses of equal amplitude $A_0 = \sqrt{P_0}$ and equal pulsewidth $\tau$, one of them a simple sine (no frequency modulation) and the other an LFM chirp, then the matched filter output $SNR$ of the chirp would be larger than that of the sine, precisely by a factor of $PCR$. I had reasoned that they both have equal energy $(E = P_0 \tau)$, but the chirp has a much narrower peak width than the sine $(\delta t = \tau/PCR)$, which would result in a larger peak height, inversely proportional to keep the total energy fixed: $$E = (P_0) (\tau) = (P_0 PCR)(\frac{\tau}{PCR}) $$

This is wrong, yet for some reason it's stated exactly this way on Wikipedia's article on the subject. Someone should fix this?

What actually happens

It turns out that both signals have the same peak power, despite the fact that the chirp's energy is squeezed into a smaller width. Therefore, we conclude that the total energy of the matched filter output of an LFM is actually smaller!

At first this looks like a problem. Moreover, you're not going to do any better with any other filter, because matched filters are optimal. But as has been noted already, the same peak power will be observed regardless of the waveform. Improved $SNR$ isn't the advantage of an LFM. The real benefit is the resolution. With the same $\tau$, I can resolve echoes that are less than $\tau$ apart, without losing $SNR$. That's pretty awesome, it turns out.

How to actually observe "Pulse Compression Gain"

If you really want to see pulse compression gain, then you've gotta compare something different. Indeed, the subject of pulse compression is usually introduced by explaining the tension between time resolution and $SNR$, before pulse compression. With a simple sine, you need to have a longer $\tau$ to improve $SNR$ (more total energy), but then you have the problem that you can't resolve echoes as well as you did before. If you want to resolve echoes better, you have to either crank up your peak power or sacrifice $SNR$.

As I mentioned above, pulse compression allows you to resolve shorter time-scales while keeping $SNR$ fixed. Alternatively, you could try comparing two signals with equal bandwidth: a simple sine wave with pulsewidth $\tau$ and an LFM chirp with pulsewidth $\alpha \tau$. Both pulses have the same resolving capability when $\alpha = PCR$ of the chirp. But the chirp will a larger peak power, and thus larger $SNR$, by a factor of precisely the $PCR$. This is simply because you have a longer pulse, not because more energy has been "squeezed" into a smaller time bin.


Reverse Engineering Bessel LPF (sample frequency needed) I have the following sourcecode. For all I know, it should be a Bessel LPF at 10Hz. Is there a way to get the sampling frequence out of this?

double LPF(register double value) {
    static double buf[4];
    register double tmp, fir, iir;
    tmp= buf[0]; memmove(buf, buf+1, 3*sizeof(double));
    iir= val * 1.865660386328251e-06;
    iir -= 0.9048452950692862*tmp; fir= tmp;
    iir -= -1.898666438223082*buf[0]; fir += buf[0]+buf[0];
    fir += iir;
    tmp= buf[1]; buf[1]= iir; val= fir;
    iir= val;
    iir -= 0.8712702520352774*tmp; fir= tmp;
    iir -= -1.866439169696869*buf[2]; fir += buf[2]+buf[2];
    fir += iir;
    buf[3]= iir; val= fir;
    return val;
}

I don't think Deve's answer is correct. The coefficient determine the transfer function in normalized frequency. You can determine at what relative frequency the transfer function hits -3dB. Let's say it's 0.01. Since you  know this is corresponds to 10 Hz you can infer that the sample rate is indeed 1000Hz.
This code is difficult to read and there isn't a single line of comment in it, nor is there any explanation of what the state variables so it's difficult exactly which number is what coefficient (b0,b1,b2, a0,a1, a2 etc) but it is in principle possible to infer the sample rate.


Selective image compression? Forgive me if this is the wrong site for this question. Here's my problem: I am creating a game where players have to search for a hidden detail within the overall image. I downloaded an 11MB compressed photo from the internet which is 85MP in size (14000 x 6000 pixels). I had a very small image within this (a pixel-sized number). The color and brightness of this number closely matches the surrounding pixels so that it blends in. Now, when I went to save this image, the compression destroyed the hidden number, even at JPEG 12/Maximum. This maximum compressed JPEG increased the file size to 50MB even though the original was just 11MB. An uncompressed save file (TIF/PNG) reaches up to 250MB.

Is there anything I can do to compress the image while still preserving the integrity of the number? The number only encompasses a few pixels, but the addition of this small data has increased the file size 25x. I have access to Photoshop and Gimp. Thanks in advance, if you need to see the file then I can provide a link to it (at 250MB size or 50MB size).

To insert a JPEG image into another without recompression, you can use a jpegtran version with crop 'n' drop. It is safest to insert into coordinates that are multiples of the 16x16 pixel JPEG block size. On command line:

jpegtran.exe -drop +16+16 pixelated.jpg piano.jpg out.jpg

pixelated.jpg

piano.jpg

out.jpg

Difference images made in Gimp show only black, so there was no lossy recompression.

Difference image: piano.jpg vs. out.jpg

Difference image: pixelated.jpg (aligned) vs. out.jpg

Before all this, you can use

jpegtran -crop

to copy the target region in the big image to a file for modification.

I wonder if jpegtran can handle your huge files. My test files were quite small.


Thresholding a difference image with a very high kurtosis I have an difference image produced by substracting two images of a certain location taken a two different time steps. The resulting image of this operation shows a histogram with a very large peak centered a 0. Since this naturally can be identified with no-change between dates I interpret this as a difference image showing a very small amount of change. Yet by the magnitudes of the observations located at the tails; strong change. I would like to threshold this image to identify three classes: negative change, no-change and positive change. I have been experimenting with classical thresholding techniques (Otsus, Kapurs, Rosins, etc) applied to each of the sides of this unimodal distribution to achieve this. I beleive I am not getting to good results because these techniques are not designed for a very spikey "distribution" with very long tails. Can anyone propose a method that can be suitable for this problem or some sort of preprocessing of the original data (difference image) which can help me out on this?

The first thing that jumps to my mind is that you want to test if your data is significantly different from zero. Determining if it is positive or negative should be trivial.

I would set a threshold based on the standard deviation of the results. Anything greater than say $\pm2\sigma$ from the mean (zero) are considered to be significant changes. The suitable value of sigma to choose on how peaked your distribution is and how many and how big changes you expect.

This sort of approach is traditionally used to look for outliers, conventionally with $\pm3\sigma$, which is quite similar to what you want.


Complex low pass filters I am searching for documentation on the behaviour of complex low pass filters, and their application in telecommunications.

Thank you

A low pass filter has a frequency response that satisfies

$$|H(\omega)|\approx 0,\quad |\omega|>\omega_c\tag{1}$$

where $\omega_c$ is the cut-off frequency. A complex low pass filter must also satisfy

$$ H(\omega)\neq H^*(-\omega)\tag{2}$$

which causes its impulse response to be complex-valued. So the frequency response of a complex low pass filter is characterized by Eqs $(1)$ and $(2)$.

In discrete-time, you have the usual convolution sum

$$y[n]=\sum_{k=-\infty}^{\infty}h[k]x[n-k]\tag{3}$$

where $y[n]$ is the output sequence, $x[n]$ is the input sequence, and $h[n]$ is the filter's impulse response. Note that all sequences in $(3)$ are generally complex-valued. If you split up (3) into its real and imaginary part you get

$$y_R[n]=\sum_{k=-\infty}^{\infty}\left(h_R[k]x_R[n-k]-h_I[k]x_I[n-k]\right)\\
y_I[n]=\sum_{k=-\infty}^{\infty}\left(h_R[k]x_I[n-k]+h_I[k]x_R[n-k]\right)\tag{4}$$

where the subscripts $_R$ and $_I$ denote the real and imaginary parts, respectively. From $(4)$ you see that you have to implement two filters, one with impulse response $h_R[n]$ and the other with impulse response $h_I[n]$, and that you have to filter the real part as well as the imaginary part of the input signal with both filters.

As for applications in digital communications, anywhere where you have in-phase and quadrature components, you can model these two signals as a complex signal, and any filtering on this complex signal requires a complex filter (or, at least, a filter with complex input and output). However, in practice you most often encounter two special cases:

The filter is actually real-valued but has a complex input and output, i.e. you simply filter the real and imaginary parts of the input signal independently with the same filter. In Eq. $(4)$ that means that $h_I[n]=0$. This is the case when you apply a low pass filter after complex demodulation.
The filter is complex-valued but the input signal is real-valued, so you filter one (real-valued) signal with two filters (the real and imaginary parts of the impulse response). In Eq. $(4)$ that means that $x_I[n]=0$. This occurs when you filter a real-valued passband signal with an analytic band pass filter before complex demodulation. (Note that here the filter is a band pass, not a low pass).

A complete complex filtering operation with complex-valued inputs and outputs, and a complex impulse response (as specified by Eq. $(4)$) is used for modelling a digital communication system using an equivalent discrete-time baseband channel. The corresponding baseband channel filter is generally complex-valued.


Overview Noise Reduction Techniques I am trying to learn signal processing and when I research about ways to reduce noise in different situations it's so easy to get lost in one solution or go down the wrong path

for example, is your noise correlated? Is it not? Is it from multiple sources? Is it not?

I need to get some intuition on common solutions to these scenarios, but I too easily get lost in the massive amount of really specific research papers

Can someone point me to an overview of these techniques and when to use them etc...? Any help would be useful at this point

The frequency domain conjugate multiplication (correlation) of the received signal with the reference signal followed by the power delay profile will provide you the overall signal to noise ratio. As explained in paper titled "SNR Estimation based on Sounding Reference Signal in LTE Uplink".

There are various noise reduction algorithm exist. Which will reduce the noise and will enhance the main expected signal component. There are adaptive filter which will converge to the expected value and depend on there performance, complexity is dependent. For example LMS, NLMS and RLS are the type of adaptive filters. Which will reduce the noise but the accuracy is directly proportional to the complexity.

By filtering we can reduce the noise, but we cant remove it totally for example AWGN noise.


Cross Power Spectrum Is it necessary that cross power spectrum (CPSD) should always be between power spectrum (PSD) of individual signals, irrespective of the method we compute it? The images are shown below.

Whenever talking about Cross Correlation of signals it is important to well define it.

If you mean Cross Correlation in most common sense (Signal Processing):

$$ \left( x \star y \right) \left[ n \right] = \sum_{i = 0}^{N - m -1 } {x} \left[ m \right] y \left[ n + m \right] $$

Then of course not.
You can have cases where the DFT of the Cross Correlation has higher absolute values then the DFT of the Auto Correlation.

  The Schwarz inequality answers your question, Yes.
  But don't use pwelch for the cross spectrum. The cross spectrum is complex valued.

The above by @Stanley Pawlukiewicz is wrong.
Simple way to see so is by defining $ y \left[ n \right] = 10 x \left[ n \right] $ which of course mean the cross correlation between those is higher than the cross correlation between $ x $ and itself.

The reason is that Cauchy Schwarz Inequality is about Inner Product Space.
So if we talking about signals as vectors then we need to use the Inner Product which is the cross correlation at $ 0 $.
If we talk about Random Process than we need to use their definition (Which normalize the energy and then everything works).

If you use Cross Spectrum while sticking to the Random Process definition and normalize the Cross Correlation by the variance of the processes than indeed by Cauchy Schwarz Inequality you're guaranteed to see what you got.


A signal $g(t)$ is passed through a squaring device and then through an LPF such that bandwidth of the LPF tends to 0 I am solving a question, where a signal g(t) is passed through a squaring device and then through an LPF such that bandwidth of the LPF tends to 0.

I understand that since for this filter the bandwidth, $ \Delta f \to 0$ , the transfer function is practically an impulse. I also understand that the width of the impulse should be $ 2 * 2 \pi \Delta f = 4 \pi \Delta f $. But I don't understand why the transfer function is expressed as $ 4 \pi \Delta f \delta(f) $, that is, why is it expressed as an impulse function of area $4 \pi \Delta f $

Here is the question and the official solution to the question:

  Solution:

I don't think that using a Dirac impulse really helps to see what's going on in this problem. Using the notation of the given solution we have

$$\mathcal{F}\left\{g^2(t)\right\}=A(\omega)\;\Longrightarrow\;A(0)=E_g$$

and

$$Y(\omega)=A(\omega)H(\omega)$$

where $H(\omega)$ is the frequency response of the ideal low pass filter.

The output signal is given by

$$\begin{align}y(t)&=\frac{1}{2\pi}\int_{-\infty}^{\infty}Y(\omega)e^{j\omega t}d\omega\\&=\frac{1}{2\pi}\int_{-\infty}^{\infty}A(\omega)H(\omega)e^{j\omega t}d\omega\\&=\frac{1}{2\pi}\int_{-2\pi\Delta f}^{2\pi\Delta f}A(\omega)e^{j\omega t}d\omega\\&\approx\frac{1}{2\pi}A(0)\cdot 4\pi\Delta f\\&=2E_g\Delta f\end{align}$$

with the approximation being justified for $\Delta f\to 0$.


Filtering $ n \times n $ Images by Separable $ m \times m $ Filters. Computation Time for Filtering Using FFT, 2D Convolution and Two 1D Convolutions Consider filtering square $n\times n$ images by square, separable $m\times m$ filters.

Give general equations for the computation time for the following approaches to this:

Filtering using FFT
2-D convolution
Two 1-D convolutions

Use capital letters to indicate implementation-dependent constants.

For particular implementations of the methods of filtering in part 1) the following computation timings were recorded:
n=64    m=3   FFT=2,498,561      2-D convolution=73,829      two 1-D convolutions=245,761
n=256   m=5   FFT=53,084,161     2-D convolution=1,179,649   two 1-D convolutions=3,932,161
n=1024  m=3   FFT=1,059,061,761  2-D convolution=18,874,369  two 1-D convolutions=62,914,561
n=1024  m=11  FFT=1,059,061,761  2-D convolution=471,859,201 two 1-D convolutions=314,572,801

Estimate the values of the implementation-dependent constants to 1)
For the implementation in 2) give values of $m$ and $n$ for which the FFT implementation is expected to be fastest.

My point of view about:

FFT is: Firstly, applying FFT to each column separately and then to each row. FFT for $N$ data points takes about $2N\log N$ multiplications, for $n\times n$ images it will be $2n\cdot 2n\cdot\log n$ multiplications. Secondly, inverse FFT to get back to sensible image. Another $2n\cdot 2n\cdot \log n$. Thirdly, multiply by FT of the kernel, it will be $n^2$ complex multiplication, that is $4\cdot n^2$ real multiplication. Thus the computation time for FFT is $\mathcal O_1\left(8\cdot n^2\log n +4\cdot n^2\right)$. Is that correct?
2-D convolution: I am confused about this. It could be $\mathcal O_2\left(n^2\cdot m^2\right)$. Is it correct?
Two 1-D convolution: Firstly, convolution per row. It should be $(n+m-1)^2$. Secondly, transpose the matrix of the image. Will this process include the computation time? Thirdly, convolution on the rows again. Another $(n+m-1)^2$. Fourthly, second transpose, return the image back to its original orientation. Computation equation is $\mathcal O_3\left(2\cdot\left(n+m-1\right)^2\right)$.

However, according to my computation equations, I cannot compromise my equations to the recorded computation timings. In this way, I think I need to correct my equations to solve the problems in part 3).

I'm not sure exactly what you're after, but just to try to add data:

Using the separability property of the filter is always the right choice.
Given it is separable we now have to apply 1D convolutions. Selecting which method to do it (Frequency Domain / Spatial Domain) depends on the length of the filter and the number of pixels in each dimension of the image.
There are cases where it is better to do FFT on the Rows and Spatial Convolution on the Columns.
This is all true when "Counting" FLOPS. Real life timing is more than that. for instance, if you're using highly tuned Convolution implementation and yet "Classic" DFT implementation you might be faster doing the Spatial way even for dimensions the analytic calculation predict it is better to be done in frequency domain.


Why isn't bilateral filtering shift-invariant? I have a source (that I can't share) that states that bilateral filtering isn't shift-invariant (that is, it depends on the input signal), hence we cannot use convolution and the Fourier convolution theorem (that is, convolution is the spatial domain corresponds to multiplication in the spectral one).

Why isn't bilateral filtering shift-invariant? Why exactly can't we apply convolution and the convolution theorem to the bilateral filtering of an image? Why wouldn't a Gaussian filter also be signal dependent?

I understood that bilateral filtering combines domain and range filtering, that is, when updating a pixel of an image, it weights the contribution of other pixels in the image using both their Euclidean distance (to the pixel being updated) and their intensity (in case of grayscale) value. However, I don't get why we couldn't also use the convolution theorem to filter an image using BF.

The bilateral filter can be written as:

$$y[i, j] = \frac{\displaystyle\sum_{k,l\in \mathbb Z^2}x[i+k, j+l]\,f\Big(x[i+k,j+l] - x[i,j]\Big)\,g\Big(k^2+l^2\Big)}{\displaystyle\sum_{k,l\in \mathbb Z^2}f\Big(x[i+k,j+l] - x[i,j]\Big)\,g\Big(k^2+l^2\Big)}\tag{1}$$

Shift invariance can be defined as any shift $m,n$ of the input $x$ resulting in the corresponding shift in the output $y$:

$$x[i,j] \to x[i-m,j-n]\quad\Rightarrow\quad y[i,j]\to y[i-m,j-n].\tag{2}$$

By doing the first substitution on all occcurrences of $x$ on the right side of Eq. 1 we get as the output:

$$\begin{gather}\frac{\displaystyle\sum_{k,l\in \mathbb Z^2}x[i+k-m, j+l-n]\,f\Big(x[i+k-m,j+l-n] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)}{\displaystyle\sum_{k,l\in \mathbb Z^2}f\Big(x[i+k-m,j+l-n] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)}\\
=\frac{\displaystyle\sum_{k,l\in \mathbb Z^2}x[i-m+k, j-n+l]\,f\Big(x[i-m+k,j-n+l] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)}{\displaystyle\sum_{k,l\in \mathbb Z^2}f\Big(x[i-m+k,j-n+l] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)},\end{gather}\tag{3}$$

where we can recognize $i-m$ and $j-n$ as $i$ and $j$ of Eq. 1, so we have:

$$= y[i-m, j-n],\tag{4}$$

showing that the bilateral filter is shift-invariant.

However, these are shifts with integer $m$ and $n$. The non-linearity of the bilateral filter will make it not invariant to fractional shifts carried out by (sinc) interpolation. Perhaps this is what your source was referring to.


Should I sample at twice the bandwidth or twice the highest frequency? I am confused about what rate I should sample at.  I've heard 2 different ways:
1) Sample at twice the highest frequency
2) Sample at twice the bandwidth

If I have a signal composed of just cosine, such as:
$$x(t) = \cos(4 \pi t) + \cos(2 \pi t) + \cos(1 \pi t)$$
What is the correct minimum sampling rate?

Well, according to Nyquist-Shannon theorem you should sample at a rate which is at least twice the highest frequency you want to capture. This is also referred to as the sampling theorem because it "forms the basis" for sampling.

Now regarding your specific question, I have to say that cos(4*pi) (along with the rest of the components) is just a number which is equivalent to a DC signal. In case you imply cos(2 * pi * 4 * n/fs) (the same for the rest of the components), where n is the sample index (constantly increasing) and fs is the sampling frequency, you should sample at twice the highest frequency (which in this case is 4), which is 8.

As you might be able to say, this would produce a result (just for the first component, but similar for the rest) cos(n * pi), which has a period of 2..., which in turn means that you would get precisely two samples per period. Now in the case you happen to sample at the zero cross, all samples would end up being zero (just for this highest frequency).

I strongly suggest you have a look at some academic literature regarding sampling. I would suggest Proakis' & Manolakis' Digital Signal Processing - Principles, Algorithms and Applications and Orfanidis' Introduction to Signal Processing


What is Hue and Saturation? I have difficulty in understanding what is Hue and Saturation. please explain what is hue and saturation to me with an example or analogy.

Hue is the main indication of color. It is the value actually telling you which color it is or the value that lets you go "red" when you see a red object.

Saturation is the perceived intensity. In other words it is a value of how dominant the color is, or how colorful the object looks.

Practical example: Regardless of the color, shadows generally happen to have low saturation values.


Histogram Counts of a 2d Matrix I am currently working on a project in which I have a 10 x 2 matrix: A. I want to find the top 3 number of occurrences of each row using MATLAB.

 A = [ 1 2; 1 2; 2 4; 3 2; 3 2; 2 4; 1 2; 1 2; 5 8; 3 2]

In return I want to have:

[ 1 2 ], 4 times.
[ 3 2 ], 3 times.
[ 2 4 ], 2 times.

I have already tried hist3 and similar functions but couldn't get what I wanted.

P.S. My original matrix of statistics is very large (around 565000 x 2 double, depending on the image I analyze) and I am looking for around top 30 most occurred rows of it. In that case, I hardly believe if any loop based answer will cover it.

Edit: I have found the answer yesterday with the help of my professor.

What I have basically done is to set the range for hist3 on x and y axises. Let's say ( -5 : 5 ).

ctrs{1} = -5 : 5;
ctrs{2} = -5 : 5;
hist_A = hist3(A, ctrs);

Then by using sort function in descending format, I have calculated the occurrences and the indices of the elements of matrix A.

[hist_A_occurrence, hist_A_index] = sort(hist_A(:),'descend');

Finally, to pick the top 3 number of occurrences:

hist_A_index = hist_A_index(1:3);
hist_A_occurrence = hist_A_occurrence(1:3);

Performance:

It was able compute on A matrix, which is 563160 x 2, for the top 60 values in 0.147038 seconds with a Core i7 4th gen processor and 16gb ram.

The method is very quick, I hope this helps to the others.

Thank you all for the answers.

Are you looking for a specific function to perform this task? Otherwise this could probably be done with a simple for loop.

Determine your unique rows using the function

  bins = unique(A, 'rows');

Then just loop through each row of A and count each instance of each bin you find in A.

This is probably good enough if you dont care about how long the program runs (which still isnt long at all if your A matrix isn't ridiculously large), but keep in mind MATLAB doesn't "like" running loops.


How to perform spectral inversion in the frequency domain to convert a low-pass filter into a high-pass filter? To convert a linear-phase FIR low-pass filter into a high-pass filter with the same cut-off frequency, we can invert the sign of the low-pass filter's impulse response $h(n)$ and then add one to the center point: $\delta(n)-h(n)$, where the $\delta(n)$ is the unit impulse function.

I am trying to verify the result in the frequency domain. If the DFT of $h(n)$ is $H(f)$, then the DFT of $\delta(n)-h(n)$ is $1-H(f)$. However, in MATLAB it seems that I cannot generate the desired high-pass filter kernel by applying the ifft() function to $1-H(f)$. The unit impulse function is added to the first point of the high-pass filter kernel instead of the center point. What is the reason for this and how can I correct my code?

MATLAB code to generate this figure:

%% Generate window-sinc filter coefficients
% CutOffFreq and TransBandWidth should be expressed as fractions of sampling rate between 0 and 0.5.
CutOffFreq = 0.1;
TransBandWidth = 0.01;

% Approxiately calculate the filter order (must be a even integer).
M = 4 / TransBandWidth;
M = round((M-2)/2)*2+2; % Round to the nearest even integer

% Generate the sinc function. Filter length = filter order + 1.
idx = 0:1:M;
Sinc = sin(2*pi*CutOffFreq*(idx-M/2)) ./ (idx-M/2);
Sinc(M/2+1) = 2 * pi * CutOffFreq;  % Solve the divided by zero problem

% Generate the Blackman window.
BlackMan = 0.42 - 0.5*cos(2*pi*idx/M) + 0.08*cos(4*pi*idx/M);

% Generate the filter coefficients (the impulse response function).
b = Sinc .* BlackMan;

% Normalized the filter coefficients to have unity gain at DC.
b = b / sum(b);

%% Spectral inversion in the frequency domain
bfft = fft(b);
b1 = ifft(1-bfft);

figure('Color','w');
plot(b,'Color','b','LineWidth',1);
hold on;
plot(b1,'Color','m','LineWidth',1);
xlim([-1,length(b)+1]);
legend({'Low-pass filter kernel','Wrong high-pass kernel'});

The flaw in your logic is assuming that the frequency response of your lowpass filter in its passband is identically $1$. You're inherently assuming this when thinking that you can form the highpass filter by simply calculating $H_{hp}(f) = 1 - H_{lp}(f)$.

However, the lowpass filter has group delay because the peak of the sinc function doesn't occur at $n = 0$. Instead, as a symmetric FIR filter, it will have a delay of $\frac{N}{2}$ samples, where $N$ is the filter order. This means that the frequency response of the filter across the passband may have unity magnitude, but it will have linear phase! Recall the time shift property of the discrete-time Fourier transform:

$$
x[n-k] \Leftrightarrow X(f)e^{-j2\pi fk}
$$

So, instead of your prototype lowpass filter having a response of $1$ in the passband, it actually has a response of $e^{-j2\pi fk}$, which imparts a phase shift that is linearly proportional to frequency. If you would like to transform this filter to a highpass structure with equivalent group delay, you would need to take this into account when synthesizing your highpass filter in the frequency domain:

$$
H_{hp}(f)= e^{-j2\pi f\frac{N}{2}} - H_{lp}(f)
$$


Understanding Active Contour results - Williams and Shah 1992 I have implemented a version of Williams and Shah's 1992 active contour algorithm, which uses the energy function:

Ej = alpha * Econtj + beta * Ecurvj + gamma * Eimagej

where j is a N8 neighbor of a contour point in question. If the energy Ej is minimized by a neighbor, then the contour point in question is replaced by that neighbor.

At this point, the results seem somewhat reasonable, but I am seeing an unexpected sharp edge across the beginning and end points. Following from the paper, I am processing the first point twice to ensure that the contour is closed. However, I have included two GIFs of output: the first one includes the double counting of the first point, and the second one does not. The first point seems to be shooting away from the object in both cases.

The first point is filled in, and the last point is surrounded by a larger circle.

Is there any obvious intuition into why this may be happening?

Thanks!

Annnd I figured it out. It seems that it was just coincidence that the first/last point sits at an especially concave point on the object. I changed the ordering of my initial contour so that the start/end point sits over a different section of the object and the first point no longer shoots out.


Difference between compressive sensing and DCT-based compression? I am working on transmitting EEG signals over wireless body area network. I have applied two different compression techniques: DCT-based and compressive sensing (CS-based) approach. I noticed that the reconstructed signal is better using the DCT-based approach. Is this reasonable? I know that compression using DCT is lossy. The metrics used for evaluating the quality of the reconstructed signal is the normalized mean square error (NMSE) and the structural similarity index metric (SSIM).

Thanks in advance

In case of not fulfilling the CS theoretical conditions, is it recommended to use DCT based compression for compressing biomedical signals to transmit them wirelessly over a body area network?
I am asking this question because most papers recommend the use of CS compression, due to the energy efficient consumption property

Comparing compressive sensing and traditional compression is a complicated task. Even  comparing two compression algorithms.

One the one hand you have compression ratios, making sure the bits sent can be decompressed without other hidden bits. On the other hand, error measurements, rarely perfectly matched to the actual meaningful distorsions. [EDIT] Standard SSIM or NMSE are not well adapted to sparse signal.

Then, knowing the range over which one algorithm is better than the other is important, because one can be better at low-rates, the other at high-rates. There is rarely a free lunch in compression.

From a very mundane point of view, if you are in a compression range where the signal is not "sparse" enough (with respect to some basis) and the CS quantization is not really taken into account, yet sufficiently harmonic, it is well understandable that a DCT-like engineered algorithm can do better than a CS method if the CS theoretical conditions are not fulfilled.


Non-orthogonal basis which offers better energy compaction in image compression I am aware of image-independent basis, i.e. DCT, and image-dependent basis, i.e. Karhunen–Loève, which are used in compacting energies for image compression. These bases are orthorgonal.

Are there any compression basis which are non-orthogonal that delivers better compaction than the optimal KL transformation? They can either be independent or dependent upon the statistics of the image.

EDIT: I have forgotten to mention that I am looking at block-based transformation coding.

The motivation of the question is due to the following paper Optimized Nonorthogonal Transforms for Image Compression. The premise of this paper is that the block based Karhunen–Loève transformation yields compacted coefficients, but the coefficients are still correlated across different blocks as it is often the case where adjacent blocks of a natural image are correlated to one another. To eliminate correlation across blocks, the optimal transformation would involve the use of non-orthogonal basis.

Given a signal $x(t)$ of energy $\mathcal E$, suppose that $x(t)$
can be expressed as $\sum_i x_i\psi_i(t)$. If the $\psi_i(t)$ are a
complete set of orthonormal signals, then
$$\sum_i |x_i|^2 =\mathcal E$$
and one measure of energy compaction is the number of $x_i$ whose
squared magnitude is larger than some fraction of $\mathcal E$,
say $0.01\mathcal E$.  Another could be the smallest number of
$x_i$ such that $\sum_i |x_i|^2$ captures $99\%$, say, of the energy.

Some sets of orthonormal signals can
be much better for representing a specific $x(t)$ than other sets of orthonormal
signals. For example, sinusoids and complex exponentials are commonly
used to represent periodic signals, and do a good job with respect to
energy compaction when used with low pass or bandpass signals.  But they
require infinitely many components to represent signals such as
sawtooth waves that can be represented using just one or two coefficients
using orthonormal Legendre polynomials.
So why are generations of electrical engineering students forced to
calculate Fourier series of sawtooth waves at all?  Well, Fourier
series simplify life
in other calculations such as what happens when the signal passes through
a linear time-invariant system while Legendre polynomials do not
If signal-dependent orthonormal bases are allowed, then if we start
a Gram-Schmidt orthonormalization procedure beginning
with taking $\psi_1(t)$ to be the unit energy signal
$x(t)/\sqrt{\mathcal E}$ as the first member of a
set of orthonormal signals, then this orthonormal set can be used to represent $x(t)$ with
$x_1 = \sqrt{\mathcal E}, x_i = 0, i > 1$ for maximal compaction!  In fact,
we don't even need to bother to find $\psi_i(t)$ for $i > 1$.
So, how much compaction can be achieved with signal-dependent
bases depends on the degree of dependence allowed between the
signal and the basis to be used, etc.
The representation of $x(t)$ with respect to a set of non-orthonormal signals
is more complicated to determine.  With orthonormal signals, the coefficient
$x_i$ of $\psi(t)$ is easily determined since
$$\langle x(t), \psi_i(t)\rangle = \langle \sum_k x_k\psi_k(t), \psi_i(t)\rangle
= \sum_k x_k \langle \psi_k(t), \psi_i(t)\rangle = x_i,$$
since
$$ \langle \psi_k(t), \psi_i(t)\rangle =
\begin{cases}1, & \text{if}~k = i,\\0,  & \text{if}~k \neq i.,\end{cases}$$
while the simplification obtained via the orthonormality is not available
for non-orthonormal signals.  The question of how many of the $x_i$ are
needed to capture most of the energy is also difficult to determine
since the energy is not $\sum_i |x_i|^2$; it is necessary to take
the cross-products $x_ix_k$ into account as well.

So, unless there are specific
reasons why a particular non-orthonormal basis is important to use
(e.g. it helps with things one needs to do in image processing
or video processing or audio processing etc), simply using such a
basis because it provides better energy compaction might well be putting the cart
before the horse.  But, as always, your mileage may vary, and there may
be non-orthonormal signal bases that not only achieve better energy
compaction but also ease your life in other aspects.  It is just that
it is doubtful that a general theory exists saying that non-orthonormal
bases achieve better energy compaction than orthonormal bases.


How to determine degree of polarization in time domain! I know there are methods to determine the degree of polarization of multi-dimensional signals, say [Bx(t), By(t)] -> DOP(f). Is there any way to extract a DOP for given vector time series in the time domain? Thanks!

You should probably use one of the algorithms for primary component (PCA) or independent component (ICA) investigation. They heavily relay on covatiance matrices and their eigenvalue-eigenvector decomposition.

If your time series represent an oscillation with dominant direction, the eigenvector for the maximum magnitude eigenvalue would represent that direction. If the signal is some kind of elliptic, the basis of eigenvectors would represent the oscillation ellipsoid dimensions. Degree of polarization is then easily yielded from relations of main-to-other eigenvectors magnitudes.


Don't understand what is meant by signal dimension I don't understand the concept of dimension of a signal.  I ran into it in an explanation of Shannon Capacity, and in a paper on spread spectrum. I was hoping somebody could explain with an example.  Does it apply to analog as well as digital?  For instance is AM two dimensional considering time and amplitude?  Or is QPSK two dimensional because it is a combination of a sin and cosine term?  Or does multilevel signaling like Pulse Amplitude Modulation have as many dimensions as discrete levels it can take, or would dimensionality refer to the number of possible symbols it can represent in a pulse?

In digital communication systems, the dimension of a modulation scheme refers to the number of basis function, i.e. independent/orthogonal signals (in this case $\sin$ and $\cos$ function) used to represent the symbols. In these systems $k=\log_2 M$ binary digits are mapped into analog waveforms below
$$\left\{s_m(t), m = 1, 2, \ldots, M \ \big\vert \ M = 2^k\right\}$$
Take a PAM signal for instance, this has the form
\begin{align}
s_m(t)&=\Re\left[A_mg(t)\exp\left(j2\pi f_c t\right)\right]\\
&= A_mg(t)\cos\left(2\pi f_c t\right)
\end{align}
One basis function (one axis), the in-phase component (here $x$-axis) is used to represent the signal. The signal is one dimensional. Now consider PSK modulation, the signal waveforms are represented as
\begin{align}
s_m(t)&=\Re\bigg[g(t)\exp\left(j2\pi\frac{m - 1}{M}\right)\exp\left(j2\pi f_c t\right)\bigg]\\
&= g(t)\cos\left(2\pi f_c t + 2\pi\frac{m - 1}{M}\right)
\end{align}
Which can be decomposed down into sine and cosine components. This is a two-dimensional signal. The same goes for QAM signals.

A common misconception is taking the dimension for the number of bits per symbol. In brief, the dimension is the number of orthogonal basis functions used in the linear combination representing the waveforms $s_m (t)$. Note that despite the digital part of it, what is not always aparent is that the modulator takes in charge the mapping of the binary info into the symbols in the analog domain ready for transmission.

Having said that, there is also a notion of multi-dimensional signals where the time and frequency intervals are subdivided into short intervals for transmission. More on this and my explanation above see $[1]$.

$[1]$, J. G. Proakis, Digital Communications, 4th edition, McGraw Hill, chapter 4.


does audio feedback consist of harmonic content? So sound from musical instruments as well as human voice consists of harmonics. just wondering if feedback from speakers also have harmonics component
I have heard that feedback does not consist of harmonics, whereas from other sources, it says that it consists only odd harmonic frequency content. so which is true?

The contradicting statements may be a result of looking at two different aspects of the situation.  As Hilmar pointed out in an earlier answer, any clipping (be it electronic or mechanical) will indeed produce harmonics.

However, feedback distortion (not necessarily reaching to the point of the loud feedback oscillation being produced) definitely does not consist of harmonics, or even any spurious frequencies added!  Feedback distortion comes from a frequency response with a very weird structure (very unnatural to the ear, often described as adding a metallic quality to the original sound).

The reason is the following: since the sound takes a fixed amount of time to travel from the speaker to the microphones (simplified analysis assuming single speaker and single microphone, no echos/reflections, and no attenuation of sound w.r.t. distance traveled), then the frequencies at which the feedback is positive, thus producing peaks in the frequency response, are the frequencies $f$ such that $$d = k \frac{\,c\,}{f}$$ where $d$ is the distance between the speaker and the microphone, $c$ is the speed of sound, and $k$ is a positive integer.  In other words, positive feedback (and thus larger gain) occurs at frequencies where the distance between the speaker and microphone is an exact multiple of the corresponding wavelength.

Now, you see, that means positive feedback occurs at frequencies $$f_k = k \frac{\,c\,}{d}~~~~k = 1, 2, 3, \cdots$$

There, you can see that the distortion is not harmonic in nature --- first of all, the above, assuming gains are low enough that you don't produce clipping or oscillation, is linear distortion (strangely enough, seeing how the feedback distortion is a very unnatural-sounding artifact); and second, the spectral locations of the "spurious" frequencies are fixed and equi-spaced, independent of the input frequency. In fact, harmonic distortion applies to an input tone (a sinusoid at a fixed frequency), whereas feedback distortion simply introduces a "weird/unnatural" frequency response, so it applies to signals with a spectrum covering ranges of frequencies.

A feedback oscillation may be subject to a combination of both effects: it may still carry a perfectly audible "feedback distortion" quality (especially if the gain is just enough to produce the oscillation), while it may or may not produce clipping that adds harmonic distortion to the mix.


How Come the Low Pass Filter in Sobel Operator Isn't Normalized? I am relatively new to the field of computer vision and I have just learnt about the sobel operator. The sobel operator in the x direction is a convolution of the finite difference kernel $[1,0,-1]$ and the gaussian smoothing kernel $[1,2,1]$. Why is it the case that the smoothing kernel does not need to be normalised ?
For example, the vector below convolved with the image will result in pixel intensities that are higher than the original values. Eg, $[50,100,50]$ will result in the middle pixel getting a value of 300 which is not the intended effect of smoothing. If normalisation is applied, then the middle pixel would get a value 75, which blurs the image.
\begin{bmatrix}
1 \\
2 \\
1
\end{bmatrix}
I hope my question was clear in the sense that i don't see how applying $[1,2,1]$ filter results in blurring without normalisation.
EDIT
How the sobel operator is obtained.

The answer is simple, the Sobel Filter is a composition of Lows Pass Filter (LPF) and High Pass Filter (HPF). The composition is done by convolution.
Now, indeed the LPF presented above $ {\left[ 1, 2, 1 \right]}^{T} $ has amplification in the DC value (Its sum is 4 so the amplification is 4). Yet it is convolved with an HPF filter which rejects the DC component.
Convolution is multiplication in the Frequency Domain, since we multiply 4 by zero we essentially get zero.
Actually multiplication of LPF and HPF gives a Band Pass Filter (BPF) (In case they have some overlap in Frequency Domain). Hence in the case above, the Sobel Filter is actually a BPF.


Question regarding AC power of ergodic process We know Ergodic process is the subset of Weakly stationary process which permits us to substitute time average for ensemble Average

My teacher said If $X(t)$ is Ergodic random process then following relationships holds true and we even solved many questions based on these relations in which either graph of ACF(autocorelation function) or expression of ACF or PSD (fourier transform of ACF) is given and we're asked to find AC power,total power,Dc power etc.

but i don't know how to derive these expressions by exploiting property of ergodicity in mean and Auto-corelation .
$\lim_{\tau \to\infty}R_{XX}(\tau)=\text{d.c power}$

and $R_{XX}(0)$ = $\text{Total power}$

so, $R_{XX}(0)-R_{X}(\infty)=$$\text{A C power}$

any help or referance material will be appreciated.thank you

Considering real valued WSS processes whose auto-correlation and auto-covariance functions are defined as
$$ r_x(\tau) = E[ x(t)x(t+\tau) ] $$

$$ c_x(\tau) = E[ (x(t)-\mu_x)(x(t+\tau)-\mu_x) ] = r_x(\tau) - \mu_x^2 $$

Then the following are just basic observations.

The total power of a WSS random process is given by
$$ P_x = E[ x^2(t) ] = E[ x(t)x(t) ] = r_x(0) $$
For a nonperiodic, WSS, RP to be ergodic in the mean, a necessary and sufficient condition is such that the auto covariance function $c_x(\tau)$ goes to zero as $\tau$ goes to $\infty$ :
$$ \lim_{\tau \to \infty}  c_x(\tau) = 0$$
Since auto correlation sequence and autocovariance sequences are related by
$$ r_x(\tau) = c_x(\tau) + \mu_x^2 $$

then the value of ACF as $\tau$ goes infinity yields the DC power as
$$ \lim_{\tau \to \infty} r_x(\tau) =\lim_{\tau \to \infty} c_x(\tau) + \mu_x^2 = 0 + \mu_x^2 $$

hence $r_x(\infty)$ is the DC power for an ergodic process.

The difference between the total power and the DC power is the AC power.

For a zero-mean RP, the total power and the AC power are the same.


Rules of image's separability We know separable convolution helps reduce number of calculations during convolution operation. As usual we see this separability is being applied on convolution kernels. Can't we apply this separability to images as well? If we can do that, therefore, it would be a mathematically valid operation and the total number of calculations will be massively reduced. But it is unfortunate that not all 2D matrices are separable to make 1Ds. So what are the rules for a matrix to be separable and is there any mathematical formulation for that to find the separated form?

A separable image $I$ that could be written in a separable way with vectors $a$ and $b$ as:

$$ I= a^Tb$$

would necessarily be of rank one at most. Either a zero-matrix (rank 0), or with only one linearly independent row or column. This are rather dull images alone. But this could be a model for a slowly varying background, some people use that for background source separation in image sequences. The following images are such images, generated the following way. $a$ and $b$ are uniform random sequences, $a_s$ is $a$ with sorted values, and  $b_s$ is the cumulative sum of vector $b$:

However, splitting an image in a weighted sum of rank-one matrices is at play in SVD decomposition. Finding low-rank approximations to matrices is a very current problem, and depends a lot on "how close is the approximation" (ie: with what norm?).

Wikipedia: Low-rank approximation
Bart Wronski: Separate your filters! Separability, SVD and low-rank approximation of 2D image processing filters


Scaling a wavelet in continuous wavelet transform The mexican hat wavelet, obtained from the second derivative of a Gaussian, has a functional form of

$$ \frac{2}{\sqrt{3 \sqrt{\pi} \sigma \ }} \left(1-\frac{t^2}{\sigma^2} \right) e^{-t^2/(2\sigma^2)} $$

where "$\sigma$ is the standard deviation of a Gaussian and $t$ is the independent variable.

If we graphically wish to graphically see its dilations, say at scale $1 \le a \le 5$, the typical form of a wavelet function $\psi$ has a scale parameter $a$ and a translation parameter $b$. The functional form of the Mexican hat given in MATLAB does not explicitly have them both.

The standard deviation $\sigma$ is the only variable here that controls with width, but this is not strictly equal to the scale $a$. How should we generate scaled versions of the Mexican hat?

Thanks.

Okay, I am not quite understanding your question.  Let's start with the definition of the Gaussian, aka the Bell Curve, in its general form.

$$  f(t) = \frac{1}{ \sigma \sqrt{2\pi}} e^{ -\frac{(t-\mu)^2}{2\sigma^2} } $$

$\mu$ is the mean, and represents where the peak occurs.

$\sigma$ is the standard deviation, and identifies where the inflection points are.

Taking the first derivative:

$$
\begin{aligned}
f'(t) & = \frac{1}{ \sigma \sqrt{2\pi}}\left[ e^{ -\frac{(t-\mu)^2}{2\sigma^2} } \left( -\frac{(t-\mu)}{\sigma^2} \right) \right] \\
 & = - \frac{1}{ \sigma^3 \sqrt{2\pi}}\left[ e^{ -\frac{(t-\mu)^2}{2\sigma^2} } \left( t-\mu \right) \right] \\
\end{aligned}
$$

Taking the second derivative.

$$
\begin{aligned}
f''(t) &= - \frac{1}{ \sigma^3 \sqrt{2\pi}}\left[  e^{ -\frac{(t-\mu)^2}{2\sigma^2} } \left( -\frac{(t-\mu)}{\sigma^2} \right) \left( t-\mu \right) +
e^{ -\frac{(t-\mu)^2}{2\sigma^2} }   \left( 1 \right)  \right] \\
f''(t) &=  \frac{1}{ \sigma^3 \sqrt{2\pi}}\left[    \frac{(t-\mu)^2}{\sigma^2} - 1 \right]  e^{ -\frac{(t-\mu)^2}{2\sigma^2} } \\
f''(t) &=  \frac{1}{ \sigma^3 \sqrt{2\pi}}\left[  \left( \frac{t-\mu}{\sigma} \right)^2 - 1 \right]  e^{ - \frac{1}{2} \left( \frac{t-\mu}{\sigma} \right)^2 } \\
\end{aligned}
$$

In your example $\mu=0$ and $\sigma=1$.  In your graph, the zero is where the peak is, and the one is the +/-1 on your horizontal scale.  Notice that the second derivative is zero at these points (inflection in the original function.)

You can translate the graph by setting $\mu$.  You can stretch the graph horizontally by setting $\sigma$.

I am not familiar with wavelets, so I don't know what this $a$ and $b$ are that you are talking about.  You can definitely to another coordinate conversion if you want:

$$ t = a \tau + b $$

Gosh, but that is awfully similar to

$$ \frac{1}{\sigma} t - \frac{\mu}{\sigma} $$

that you already got.

Based on the comments, to retrofit the function call to MATLAB (or wherever):

$$ \mu = b $$

$$ \sigma = 0.63628 a $$


Calculation of the salience function (Klapuri 2006) I am reading this paper here and I wanted to check my understanding regarding the salience function:
$$
\hat{s}(\tau)=\sum_{m=1}^M g(\tau, m) \max_{k\in \kappa_{\tau, m}}\lvert Y(k)\rvert
$$

Where the set $\kappa_{\tau, m}$ defines a range of frequency bins in the vicinity of the $m^{\rm th}$ overtone partial of the F0 candidate.

We want to find the maximum of $Y(k)$ for those k. But k is a frequency. And Y is, according to Klapuri the discrete STFT of the signal. I don't understand what the STFT on a given frequency is.

The STFT represents a 2D vector so $Y(k)$ must be a vector, right?
How do we find the max of list of vectors?
Is something in my understanding wrong?

So let's say we are testing the salience of the period 0.002. We start with $m=1$, and get $g(\tau, m)$. Then what?

You are correct that STFT is a function of 2 variables (usually x-axis is time and y-axis is frequency). The author is a bit careless in calling $Y(k)$ the STFT. What he really means is that $Y(k)$ is a fixed column (i.e. fixed time slice) of the STFT. The salience function is more carefully defined by the same author in this paper titled Multipitch analysis of polyphonic music and speech
signals using an auditory mode Eq. 13, where he uses a subscript $t$ to denote a fixed time slice of the STFT.


Different amplitude in time-domain vs frequency-domain Reading the book Understanding Digital Signal Processing I am confused with the following picture:

First figure is the input signal ($x_1(n) = \sin(2\pi f_0nt_s)$, $f_0$ frequency, $t_s$ being the sampling period) at 1Hz. Second is the input signal squared ($y_1(n) = \sin(2\pi f_0nt_s) \cdot \sin(2\pi f_0nt_s)$). Third figure is the second signal displayed in frequency domain. How can it be that the amplitude is $-0.5$ when the second figure is strictly positive and the amplitude is $1$?

Note: I am a complete beginner at this.

How can it be that the amplitude is −0.5

Primarily sloppy notation in the book.

The frequency domain representation of a signal uses complex numbers, i.e. you would need two real numbers to describe it. Either Amplitude and Phase or Real and Imaginary Part. The correct way to describe this would

Amplitude is $0.5$ and the phase is $\pi$, or
Real part is $-0.5$ and the Imaginary part is $0$

The whole notion of negative/positive doesn't make sense for a complex number, so you can't think about this that way.


Time delay estimate and sampling frequency This is an extention to my other question. My sampling frequency is only about 2x my tone frequency. I want to estimate time delay of arrival using this signal. I am generating the tone signal so I know its properties. Can I use compressed sensing to reconstruct my tone wave received by the sensors and get a good time delay estimate? An answer with some matlab code would be brilliant.

A colleague of mine has investigated time delay estimation using compressed sensing in this paper: http://vbn.aau.dk/en/publications/compressive-parameter-estimation-for-sparse-translationinvariant-signals-using-polar-interpolation(56d35506-d4ae-4296-9a36-7ec95b1160c7).html
You can find the accompanying Matlab code here: http://www.sparsesampling.com/cpe/


Using a Kalman filter to predict samples between measurements (cross posted from stack overflow)
I have a dataset where measurements were taken at 10 Hz, and I am trying to use a Kalman filter to add predicted samples in between the measurements, so that my output is at 100 Hz. I have it working ok when the velocity is linear, but when the direction changes, the filter takes a while to catch up. I am new to Kalman models, so am very likely making some mistakes in my settings. See image

for an example, the red is measured data, with stepping in between measurements. The blue is the Kalman corrected.

std::vector measurements is a dummy data array I am testing with.

The main Kalman code is based on this:

https://github.com/hmartiro/kalman-cpp/blob/master/kalman.cpp

Am I correct in feeding the prediction back into the filter for every loop? Or is this line:

yPos << kf.state().transpose(); wrong?

my code is:

int main(int argc, char* argv[]) {

  int n = 3; // Number of states
  int m = 1; // Number of measurements

  double dt = 1.0/30; // Time step

  Eigen::MatrixXd matA(n, n); // System dynamics matrix
  Eigen::MatrixXd matC(m, n); // Output matrix
  Eigen::MatrixXd matQ(n, n); // Process noise covariance
  Eigen::MatrixXd matR(m, m); // Measurement noise covariance
  Eigen::MatrixXd matP(n, n); // Estimate error covariance

  // Discrete motion, measuring position only
  matA << 1, dt, 0, 0, 1, dt, 0, 0, 1;
  matC << 1, 0, 0;

  // Reasonable covariance matrices
  matQ << 0.001, 0.001, .0, 0.001, 0.001, .0, .0, .0, .0;
  matR << 0.03;
  matP << .1, .1, .1, .1, 10000, 10, .1, 10, 100;

  // Construct the filter
  KalmanFilter kf(dt,matA, matC, matQ, matR, matP);

  // List of noisy position measurements (yPos)
  std::vector<double> measurements = {
     10,11,13,13.5,14,15.2,15.6,16,18,22,20,21,19,18,17,16,17.5,19,21,22,23,25,26,25,24,21,20,18,16
  };

  // Best guess of initial states
  Eigen::VectorXd x0(n);
  x0 << measurements[0], 0, 0;
  kf.init(dt,x0);

  // Feed measurements into filter, output estimated states
  double t = 0;
  Eigen::VectorXd y(m);

  for(int i = 0; i < measurements.size(); i++) { //ACTUAL MEASURED SAMPLE

      yPos << measurements[i];

      kf.update(yPos);

      for (int ji = 0; ji < 10; ji++)  // TEN PREDICTED SAMPLES
      {
          t += dt;

         kf.update(yPos);

          yPos << kf.state().transpose(); //USE PREDICTION AS NEW SAMPLE

      }
  }

  return 0;
}

Thank you.

I'm sorry that it was not very clear to me what you were asking. But I think I have developed some sense on it as follows:

You are physically capturing 10 noisy sensor measurements per second, which I call it as the physical measurements. But for some reasons, you want to fill in between those true-physical data samples with new synthetic data, generated mathematically at a rate of 100 samples per second (10 times more dense than available with the physical sensor measurement)

The solution to this problem has different names, based on your interpretation and which method you choose to follow, such as a sampler, a  predictor, an estimator, a curve fitting or an interpolation...

It seems, you want to use a Kalman filter (an estimator actually) to generate those new synthetic samples in between the true physical measurements. And your question is whether you should feed those synthetically generated new sampes back into the Kalman filter (for residual calculations) as if the data were generated from a true physical sensor, so that the filter can use it in residual calculations to get rid of the noise present in it and refine the state estimates?

That's what I understand from your question set up and my answer to it is no you cannot. Those new data cannot replace true phsyical measurements. They do not involve new information that's called as the innovations process of a random process (that you're trying to estimate), hence they cannot be used to make refined computations or decisions alone.

In fact, the Kalman filter will only propagate the current state $X(t)$ into the next phsyical measurement point $X(t+T_s)$ by using the system dynamic matrix alone. (i.e., it will solve the differential equation numerically based on the initial condition provided from the previous state estimate, and sample this solution at the requested high data rate)


Is video considered a Tensor? A sampled single dimensional signal is a vector.
A sampled image is a matrix.
Is video considered a rank-3 tensor since it is a set of images stacked on top of each other?

yes - a video signal would be considered 3rd order tensor.

scalars are 0-order tensors, vectors are 1st order tensors, matrices are 2nd order tensor,  data volumes are 3rd order tensors and so on for multi-dimensional arrays.


Training own classifiers using trainCascadeObjectDetector for face detection I want to train a cascade classifier to face detect. I read this page. But I wont use trainingImageLabeler app, because I already have face and nonface database sized 24x24. How I use trainCascadeObjectDetector function with my own database?

Is there any example except this you know?

If you want to train the classifiers with your own database, you will only need 'trainCascadeObjectDetector' function and feed your images into the proper arguments (Positive, negative images). The output classifier will be in your 'outputXMLFilename' as in  traincascadeobjectdetector

trainingImageLabeler is helpful function for classifying positive/negative image. Try this function later when you work with larger database.


measuring how much signal corrupted by passing from a medium I am looking for a measure to report how much a signal has corrupted by passing a medium.

Suppose we have a signal 2Hz and amplitude 1. I pass it from a medium. A receiver got the signal with 2Hz, amplitude 0.5 and a few more signals.
Convolution gives me an array that say with a specific delay overlapping between input and output is maximum.
Do we have a measure to report a single scalar how much a signal has corrupted?

The most typical measure of a signal corruption (or degradation) would be an error squarred integral (or sum for discrete case), preferably normalized such as:

$$ \epsilon = \frac {\int (x(t)-x_r(t))^2 dt}{\int (x(t))^2 dt} $$

Where $x_r(t)$ is identified as the received and corrupted signal. The interval of integration choosen to fit into signal duration of interest.

This is an unweighted (or equiweighted) calculation. If you think that certain features of the error is more impotrant, more relevant, more effective than the others, then you could reflect this by modifications such as simple linear weighting or onlinear processing of the error accordingly, which would produce more meaningful results.


What is the forward error correction in 802.11a? I am trying to put together a baseband-equivalent simulation of 802.11a signals, which according to Wikipedia has a variable encoding (PSK, QAM, etc.) that have differing coding rates (1/2, 2/3, 3/4) depending on the choice of RATE bits.   I can't find any details on the forward error corrector used in 802.11a.  Could someone point me the the algorithm used in 802.11a?

From the official standard, IEEE Std 802.11a-1999, Part 11, section 17.3.5.5 Convolutional encoder, see the rate 1/2 code:

The resulting two streams A and B are interleaved bit-alternatingly.

The other rates (2/3 and 3/4) are achieved through puncturing (as defined on the page right after).

Since this code can't deal with anything that deserves the name "burst error", which in the context of OFDM means "wideband disturbance" (in this case, "interference spanning more than 3 subcarriers" is impossible to correct reliably), there's an interleaver following this. Don't forget that in your implementation.

By modern standards, this is not an overly impressive code, to be honest.


Determine fundamental period from zero crossing signal I have got a signal which consists of zero crossings over discrete time and I would like to estimate the fundamental frequency (period) from this signal in order to remove noisy samples.

The signal is about 500 to 1500 samples long and has about 10-50 zero crossings, e.g. x[93] = 0; x[183] = 0; x[244]; x[282]; x[310]; x[439]; x[502]; x[515]; x[570]; x[590]; x[640]; x[635]; x[650]; x[710]; x[740]; x[835]; x[850]; x[905]; x[915]; x[980]; x[1050]; x[1110];

The output should be zero crossings again or their fundamental period, but without the "noise". I'm only interested in the position of the "new" zero crossings.

I'm also not quite sure, how to model a zero-crossing signal for further processing.

Thanks for any ideas.

In general, you can not determine period of a signal from just the zero crossings.  If the signal has harmonic content, then a pitch detection/estimation algorithm, such as autocorrelation, might be one solution.  If the signal is more spectrally pure (e.g. close to sinusoidal fundamental, little overtone or harmonic energy), then a suitable DSP band-pass filter applied before locating the zero-crossings may help get rid of "noisy" samples.


Question based on scaling the FFT for obtaining PSD Following is the code given in MATLAB's site to estimate PSD using FFT.

Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t) + randn(size(t));

N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;

plot(freq,10*log10(psdx))
grid on
title('Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)'

I know I might lack a proper understanding of the fundamentals, but can anyone explain in order to extract the first half of xdft why does the index run from $1$ to $N/2 + 1$ and not $1$ to $N/2$? (I suppose the Nyquist frequency lies at $i = N/2$, am I right?)
Some sources mention that the square of the magnitude should be scaled by $\frac{1}{N}$, while here it is scaled by $\frac{1}{Fs*N}$. I am unable to figure out which of these two should be used.
I understand that the DC value should be left untouched when converting to single-sided spectrum. But why is the scaling by $2$ performed only till end-1 and not till end?

The exact location (index) of the Nyquist frequency in your frequency vector depends on the length of the input signal. If the length is odd the Nyquist frequency is not included; if the length is even it is included. Consequently, and assuming the signal is real, either the first $\frac{N}{2}$ or $ \frac{N}{2} + 1$ samples are needed to fully describe your signal in the frequency domain. Note that, like the first sample describing the DC value, the sample describing the Nyquist value only occurs once; it should therefore not be multiplied by two when calculating the PSD of a real signal.

The proper scaling value to be used depends on the exact definition of the DFT. I have not checked MATLAB’s definition but if you want to make sure that the used scaling factor is correct use Parseval’s theorem. Check if the signal power in the time and frequency domain are identical; if not adjust the scaling factor accordingly.


adding u(t) and r(t) to create a signal I previously asked this question in the maths community not realising that this stack exchange community existed. I have searched online extensively and haven't been able to find an explanation to how this works and would really appreciate some input.

I'm studying this signals and systems text book and I can't figure out the logic behind how $u(t)$ and $r(t)$ can be added together to form some signals.

$$ u(t) \triangleq
\begin{cases}
0, & \text{if }t<0 \\
1, & \text{if }t\ge0
\end{cases} $$

$$ $$

$$ \begin{align}
r(t) &  \triangleq
\begin{cases}
0, & \text{if }t<0 \\
t, & \text{if }t\ge0
\end{cases} \\
& = t \ u(t) \\
\end{align}$$

$$ $$

Basically, an example in the book gives this signal:

So the book says that it can be described by this equation:
$$
x(t) = u(t)-r(t-1)+2r(t-2)-r(t-3)+u(t-4)-2u(t-5)
$$

It states that the sum of the ramps $-r(t-1)+2r(t-2)$ goes upward with the slope equal to one starting at $t=2$. I can't understand how this works geometrically and haven't been able to find anything that has helped. Does anyone know how this type of modelling is done?

I'm not sure I understand the delta-equals notation, but I'll assume it means definition.

It's important to understand that you're doing a variable substitution here.  In the strictest sense, the $t$ in $u(t-1)$ is not the same $t$ as in $u(t)$.  It may be easier for you to understand if you write your definitions as functions of $z$ instead of $t$:

$$ u(z) \triangleq
\begin{cases}
0, & \text{if }z<0 \\
1, & \text{if }z\ge0
\end{cases} $$

$$ \begin{align}
r(z) &  \triangleq
\begin{cases}
0, & \text{if }z<0 \\
z, & \text{if }z\ge0
\end{cases} \\
& = z \ u(z) \\
\end{align}$$

So, $z=t$ in $u(t)$ and $z=t-1$ in $r(t-1)$, etc.  You're changing a function of $z$ into a function of $t$.

The two functions you're dealing with have a breakpoint in them where the function changes behavior.  In both cases it happens when $z=0$.  So, looking at $u(z)$ for example, $u(z)=0$ for all values of $z$ less than zero and $u(z)=1$ for $z=0$ and all values greater than zero-- just saying what the function says, but in words.

When I use $u(t)$, I am substituting $z=t$ so $u(t)=0$ for all times before $t=0$, and $u(t)=1$ from time $t=0$ onward.

Now think about what's happening in the function argument with all the $(t-1)$, $(t-2)$ values...  If you were to look at $u(t-4)$, then, it has the effect of shifting the entire function to the right in time.  In this case:
$$ u(t-4) =
\begin{cases}
0, & \text{if }t-4<0 \\
1, & \text{if }t-4\ge0
\end{cases} $$
or, rearranging the conditions a bit:
$$ u(t-4) =
\begin{cases}
0, & \text{if }t<4 \\
1, & \text{if }t\ge4
\end{cases} $$
so you have a step function that is zero for all times before $t=4$, and one for times $t=4$ onward.

With the $r(z)$ function, with $z=t-2$ you substitute for $z$ in all places, so you have:
$$ \begin{align}
r(t-2) &  =
\begin{cases}
0, & \text{if }t<2 \\
t-2, & \text{if }t\ge 2
\end{cases} \\
& = (t-2) \cdot u(t-2) \\
\end{align}$$
which is zero for all times before $t=2$ and then a line starting at $(t=2,x(t)=0)$ with a slope of positive 1 and a y-intercept of -2.

The coefficients have the effect of changing the slope and corresponding y-intercept but the t-intercept remains the same.  $2\cdot r(t-2)$ has a slope of 2 and a y-intercept of -4 but still crosses the $t$ axis at $t=2$.

The x(t) function is a linear combination of such functions, so you must add up the values at each point.  The easiest way to see this on graph paper would be to draw each term separately one above the other and then add them vertically.

With the two components you you've highlighted, $−r(t−1)+2\cdot r(t−2)$:  There is a slope of zero before $t=1$ because $u(t)$ is zero slope.  The first term starts a ramp of slope of -1 starting at $t=1$ and the second term adds to that a ramp of slope +2 starting at $t=2$.  Where those two terms overlap (everywhere from $t=2$ and to the right) the combined slope of those terms is +1.


Backward finite difference differentiation filter frequency response As title says what would be frequency response of backward finite difference differential filter, or what would be error of this differential filter, analyzed upon frequency of a signal?

In discrete-time signal processing terminology a backward (one sample) difference filter is
$$y[n] = x[n] - x[n-1]$$
which is a causal, LTI system with an FIR impulse response of
$$h[n] = \delta[n] -\delta[n-1]$$
The associated frequency response is then the DTFT of the impulse response as
$$\begin{align}
H(e^{j\omega}) = \mathscr{F}\{ h[n] \} &= 1 - e^{-j\omega} \\
&= e^{-j\omega/2}(e^{-j\omega/2} - e^{-j\omega/2}) \\
&= 2j e^{-j\omega/2} \ \tfrac1{2j}(e^{-j\omega/2} - e^{-j\omega/2}) \\
&= e^{-j\omega/2} 2j  \ \sin(\omega/2) \\
&\approx e^{-j\omega/2} j  \ \omega \qquad \text{for } |\omega| \ll \pi \\
\end{align}$$
So you have a delay of $\frac12$ sample, an additional phase shift of +90° (due to the $j$ factor), and a gain proportional to frequency (at least for low frequencies, relative to Nyquist).  Other than the half-sample delay, this is what you expect from a differentiator.
If you happen to use an $N$-sample backward difference such as
$$y[n] = x[n] - x[n-N]$$
then the associated frequency response will be
$$H_N(e^{j\omega}) = 1 - e^{-jN\omega} $$
This is similar to the first orde case, but the delay is $\frac{N}2$ samples, and the gain is $N \omega / 2$ .


Detect the beginning of an increasing signal After denoising and cleaning, I get amplitude signals like this (y-axis: dB):

On bottom left of each of these 3 graphs, you can see a noise floor (nearly "horizontal line"). This noise floor fluctuates in a range < +-1 dB.

How to detect when the real signal begins, i.e. when the signal starts increasing?

Here is what I've tried:

Let $t_i$ be the first time the signal goes higher than $noisefloor + i \  dB$. For example $t_4$ is the first time the signal goes 4 dB more than the noise floor.

Let's do a quadratic or cubic interpolation $f$ of the curve going by $(t_2, floor+2)$, $(t_4, floor+4)$, $(t_6, floor+6)$, etc.

Than let's check when the polynomial $f$ crosses the horizontal line $y = floor$.

It sometimes works, and better than all I've tried so far. But it's not perfect:

Here the place where the green interpolation crosses $y = floor$ (dashed line) for the first time is more or less what I'm looking for, but there are 2 crossings, how to select the best? Sometimes it's the first to choose, sometimes it's the second..., it's quite random.

Here the place where the green interpolation crosses $y = floor$ (dashed line) does not exist! so it doesn't work at all:

How to do better than this?

Note: a simple threshold-detection of when the signal goes higher than $noisefloor + 0.5 \  dB$ as "the time when the signal starts" is not precise enough for my application.

It is not entirely clear what sort of signal we are dealing with here, apart from the use of the "audio" tag. If the signal had a wider bandwidth then this would be closer to onset detection. But this is not what we are dealing with here. What we are dealing with here is a slow varying waveform that is considered "on", when it has emerged from some kind of background activity. This viewpoint would make it closer to Anomaly Detection or "Outlier Detection" type of problems.

There are a couple of approaches you can take here. One is to fit a model that includes an "activation" parameter and then try to see when does that parameter transits to activation and take that point as the beginning of the onset. If you set off down the path of model fitting, eventually you are going to have to train your model on a number of these curves so that it learns all the different possible "avenues to activation" there might be. For example, the third trace shows background activity, transiting to an intermediate plateau, transiting to full on activation and even in the "activated" region its slope might show further variation.

So, before you start looking at those techniques, maybe you can try the plain old technique of detecting outliers through the statistics of the signal.

As a human being (?) you seem to have a good idea of when you would like to consider this curve as being "on". Therefore, collect all samples of your signal within the "background" region and use something like a boxplot or fit a distribution to this data. The simplest example of that distribution would be a Gaussian that has a mean and standard deviation. This models your "normality" region. Any value that could have emerged from that distribution is dubious as to whether it belongs to the "background" or "activated" segments. But that is not true for all values because soon enough (as time evolves towards the right), the curve will start pushing towards extremal points of the distribution where the probability of generating such a value becomes smaller and smaller.

Putting a hard threshold there, would give you an estimate of where the "activation" region starts.

Hope this helps.

EDIT:

After more information about the problem was shared, I am more inclined to suggest one of the onset detection techniques which will work directly on the audio data.

In any case, the following (cave) illustrations might help a bit more with the earlier suggestions in this post.

The "Human being" comment comes into play when determining the point of earliest transition from "background" to "key on" (rather than doing it automatically). You use the data in the "background" part of the waveform to estimate the statistics of what it means for a sample to be coming from the "background" part and use that to determine a threshold beyond which the samples are now more likely to belong to the "key on" part.

Alternatively:

Combine many takes at similar settings by aligning them on the "Key on" slope and summarise all of this data with a series of boxplots, each one telling you the sort of value limits you can expect at each time instance. Use that information then to choose the point in time when there is deviation from the background.

(Inset images of boxplot and distribution from this and this wikipedia articles respectively.)


Circular convolution of real signal N by N/2 sized complex IFFT There are two related questions with good answers;
for calculating the forward FFT in N/2 values, and
calculating inverse FFT from N real values, now I'm in search for the missing link of packing the N/2+1 element complex signal X back to the N/2 element complex signal that it was derived in the first place for memory and arithmetic efficient inverse FFT.

Having calculated N/2+1 valued complex FFT from N real values (code from the first link):

x = [ 1 2 3 4 5 6 7 8 ]; n = 8; n2 = 4;
Z = fft(x(1:2:end) + i * x(2:2:end);
Ze = .5*( Z + conj([Z(1),Z(n2:-1:2)]));
Zo = -.5*j*( Z - conj([Z(1),Z(n2:-1:2)]));
X = [Ze,Ze(1)] + exp(-j*2*pi/n*(0:n2)).*[Zo,Zo(1)];

I have actually packed the real valued element X(n2+1) to imaginary part of X(1), saving some memory. Then I've point wise multiplied X with G representing the FFT of a same sized real signal g.

In order to finish the (circular) convolution of conv(x, g), I'd need to calculate the packed N/2 sized complex signal Z from (X.*G) for IFFT.

The first coefficient can be reverted as Z(1) = real(X(1)+X(n2+1))/2 + i*(imag(X(1)-X(n2+1)))/2;

As the other values X(2:n2-1) are originally produced from a bufferfly operation of
x(k) +- exp(j*k/N) * x(n2-n), I would hope that all the elements would be recoverable. The question is how.

EDIT

I need to recheck the formulas, but since for a pair A=Z(k+1), B=Z(n2+1-k):

                                         [ 1-s 1+s   c     c ] [Real(A)]
 [Real(A') Imag(A') Real(B') Imag(B')] = [ 1+s 1-s  -c    -c ]*[Imag(A)]
                                         [ -c   c    1-s -1-s] [Real(B)]
                                         [ -c   c   -1-s  1-s] [Imag(B)]

with invertible matrix, there really is a solution.

This is the general solution to recover N real samples using an N/2 Length IFFT. N/2+1 samples of the FFT are required for an even sequence N.

First generating a sample sequence of the N/2+1 samples of an FFT of a real signal:

x = [1 2 3 4 5 6 7 8];
N = length(x);
N2 = N/2;
fx = fft(x);
X = fx(1:N2+1);             % starting result: positive spectrum of real sequence

So above is the same result as can be alternatively computed with an N/2 FFT as already demonstrated.

Given such an result, below shows how to get the original sequence x using an N/2 IFFT. The code is self-explantory:

W = @(k,N) exp(-j*2*pi/N*[0:k-1]);                   % FFT Twiddle factors
Xe =  (X(1:N2) + conj(X(N2+1:-1:2)))/2;              % FFT of even samples
Xo = (X(1:N2) - conj(X(N2+1:-1:2)))./(2*W(N2,N));    % FFT of odd samples

newX = Xe + j*Xo;

xr = ifft(newX);                      % N/2 IFFT
newxo =  real(xr)                     % recovered odd x samples
newxe =  imag(xr)                     % recovered even x samples

Since each FFT in your formula for cross correlation is the FFT of a real signal, the product will also be the FFT of a real signal, hence the technique outlined above can be used.


The error between the theoretical value and the simulation value of imaginary part of the FFT results As we know, a complex analytic signal can be obtained through Fourier transformation of a real value function, and the imaginary part can be reconstructed by Hilbert transform of real part. But I have simulated a real value function using Matlab and performance FFT on it, the imaginary part obtained through FFT is very different from that obtained through Hilbert transformation of real part, Does anyone know why.

Judging from your question, I think that you're confusing two related but different Hilbert relations. The first one holds for causal signals, for which the real and imaginary parts of their Fourier transforms are related via the Hilbert transform. The other one holds for analytic signals (whose spectrum is zero for negative frequencies), for which the real and imaginary parts of the time domain signal are related via the Hilbert transform.

So if you have a real valued function, as described in your question, and you take its Fourier transform, then the imaginary part of the Fourier transform can only be obtained from its real part via the Hilbert transform, if the signal is causal. Also note that the corresponding time domain signal is generally not analytic (because in that case its spectrum must be zero for negative frequencies and, consequently, the time domain signal must be complex-valued).


How to modulate PSK63 with liquiddsp? I am trying to implement the BPSK-63 modulator using liuid dsp library.
Here is my code:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <getopt.h>
#include <math.h>
#include <complex.h>
#include <liquid/liquid.h>

#define SAMPLERATE 11025

int main(){
    modem mod   = modem_create(LIQUID_MODEM_BPSK);

    nco_crcf q = nco_crcf_create(LIQUID_NCO);
    nco_crcf_set_phase(q,0);
    nco_crcf_set_frequency(q, 2*M_PI*(2000.0/(float)SAMPLERATE));

    const char *message = "1010111100101111011001010100101110100"; //RYry

    std::complex<float> x=_Complex_I*0;
    unsigned int symbol=0;
    unsigned int outs=0;
    int ctr=0;
    float ous[1];
    int offset=0;

    while(1){
        if(ctr++ == (SAMPLERATE/63)) {
            ctr=0;
            char sym = message[offset]-'0';
            offset++;
            if(offset >= strlen(message))
                offset = 0;
            symbol = sym & 0x01;

        }

        modem_modulate(mod, symbol, &x);
        nco_crcf_mix_up(q, x, &x);
        nco_crcf_step(q);

        ous[0] = x.real();
        fwrite((void *)ous,1,sizeof(float),stdout); // | aplay -f FLOAT_LE -r 11025 -c 1
    }
}

It looks like it works, however it sounds a bit different (has some clicks) to working modulator from FLDIGI program.
However when I feed the audio into the fldigi, I can see a carrier in the right spot, however it demodulates a series of "ot e" like:
ot e
ot e
ot e

instead of RYryRYryRYryRYry.
So, it sees five characters instead of four and detects a space and a newline when there is none.
Also, i see some harmonics, but i don't think this is a problem, since fldigi provides an adjustable filter on it's waterfall.
What i'm missing?

Posting a standard is actually a good point.
After reading several definitions i've figured out that i have to keep the symbol (and thus the phase) if I want to transmit 1 and reverse it if i want to transmit zero.
(see: http://www.arrl.org/psk31-spec, PSK-63 has just higher speed, about 63 baud.)
if(ctr==(11025/63)){
    ctr=0;
    c=message[subsym]-'0';
    subsym++;
    if(subsym>=strlen(message)){
        subsym=0;
    }
    if(!c) // here instead of just assigning bit=symbol
    sym_in = !sym_in;
}


MATLAB, Removing Certain Elements From An Array That Fits The Condition removing elements from an array that have difference of 3 or less between the next element...
eg. A=[3 6 10 14 17 20] to this B=[3 10 14 20]...I believe I can use diff but in what context should I use it to achieve this ?

This code removes all elements that have diff 3 or less to the previous element.

in = [3 6 10 14 17 20 24 25 30 34 38]
in1 = in(2:end);
D = diff(in)
out = [in(1) in1(D > 3)]

output:

in =

     3     6    10    14    17    20    24    25    30    34    38

D =

     3     4     4     3     3     4     1     5     4     4

out =

     3    10    14    24    30    34    38


Windowing before viewing a filter frequency response? It is well known that it is a good idea to apply windowing to a signal, before obtaining its DFT and viewing its frequency response.

But is it a good idea to apply windowing to the impulse response of a filter, before viewing its frequency response?
Does the answer to this question depend on whether the filter is FIR or IIR?

I could perhaps see why windowing could be better than just truncating an IIR filter. But I can't decide whether the same applies to an FIR filter.

However, even in the case of an IIR filter, we usually don't do this. Instead we usually just plot the response as a function of the IIR filter's transfer function numerator & denominator coefficients (e.g. using freqz in Python or Matlab).

So is it ever a good idea to window the impulse response of a filter?

No, if you know the filter coefficients it's never a good idea to use a window. Why should it be? Instead of plotting the actual frequency response of the filter, you would plot the convolution of the filter's frequency response with the frequency response of the window, which would result in a "smeared out" frequency response. Especially for sharp frequency-selective filters the result would be misleading.

Note that the frequency response of an IIR filter is not computed by applying an FFT to its truncated impulse response (because that would be windowing with a rectangular window), but by computing the ratio of the FFTs of the numerator and the denominator. That's what Matlab's freqz does.

However, if you don't have the filter coefficients of an IIR filter, but just its impulse response, then you don't have much choice, and you have to use the impulse response. You must include as many terms as can be distinguished from noise and use an FFT. This corresponds to windowing with a (very long) rectangular window. I don't see much advantage of using any other window because you have to include practically all terms you can distinguish from noise in order to keep the truncation error small. See also the discussion in the comments below.


Why Wavelet based Transform Is More Suitable for Image Compression Compared to DCT? Everyone knows that JPEG2000 can obtain much better compression result than JPEG. The fundamental transform behind JPEG2000 is wavelet while the basic transform in JPEG is DCT. So my question is why wavelet is more suitable for image compression compared to DCT. Thanks.

Both JPEG and JPEG 2000 use the change of basis compression type.

Namely, we transform the data into a different representation assuming in this representation the number of parameters needed to describe to data is lower.
Or to the least, most of the information is gathered within few parameters.

Now, if you look at the energy level of the DCT coefficients of real world images you'd see most of the energy is limited to the very few coefficients.
This is what JPG does, keeps only the few dominant coefficients and throws the rest.

How efficient is DCT?
It depends on how less coefficients are needed to describe the image in an acceptable quality.

It turns out that for Wavelets the situation is better.
Namely on real world images less coefficients are needed to describe the image with the same perceived quality.

This property is called the ability to decorrelate the data / being energy densed, etc...

At the end, it is all up to how much less parameters are needed to describe the same data.

P.S.
Using SVD / KLT (They are equivalent in this sense) would be even more efficient in the sense of how much less coefficients are needed for the same quality (In the $ {\ell}_{2} $ sense).
The problem is those bases are adaptive to the data and hence if you use them you have to include them to be able to decompress the data.
What's so nice, you can ask which, in general, approximate the SVD better?
It turns out that Wavelets are better doing so.


What's the difference between Channel shortening and time-domain equalizer I'm asking about the difference about channel shortening and time-domain equalization. I have checked on-line and tried to get something clear but I couldn't.

What I found is that channel shortening means almost  "Calculating the coefficients for DFE equalizer" But I couldn't also imagine that.

Could you please explain.

Thank you

Yes,I agreed these terms are very common in communication engineering,

I was asking about channel shortening long time ago, and I got an answer telling me that channel shortening is just calculating the coefficients of the DFE. you can read the details in the below link

Channel Shortening for underwater acoustic channel communication

thanks


What is mean by [1 -1] Filter? While studying about Wavelet Transform, I'm confused as why this is called a [1 -1] Filter , Here is what i'm talking about.

  In this hypothetical example the student does fairly well the first
  half of the     term then neglects his or her studies for the last
  half. Thus the exam scores     for the term were 80%, 80%, 80%, 80%,
  0%, 0%, 0%, and 0%*     We can tell the average of all the scores
  (40%) and when the scores “tanked”     after the 4th exam just by
  looking. Knowing the answer in advance, however,     is a good way to
  learn and to verify the wavelet transforms. Then we     can use them
  with confidence on real-world data where we can’t simply     “eyeball”
  the final values.

  We will now walk through the CWT process step by step using the
  simplest     of the wavelet filters on this example. We begin by
  comparing the humble Haar wavelet filter, [1 –1], with the data as
  shown

80 80 80 80 0 0 0 0

If We now keep subtracting each value from the next value we will eventually get,

[0, 0, 0, 80, 0, 0, 0]

Now If I want to stretch the Filter to three points i.e If the filter is stretched from [1 –1] to [1 0 –1].

How will I calculate these points as I did in case of two point Filter where i simply subtract each value from the next value ?

I have the same book as you. :-)

The [1 -1] filter is a a simple differentiator, also known as, the Haar mother wavelet. You need to understand convolution to understand what he is saying. On the first level, you convolve your signal with [1 -1]. Then on the next level, you convole your signal with [1 0 -1], etc. That is all the CWT is doing. (Look at convolution first). If you understand that well, you can understand CWT very easily.


What does 'z' in Z-transform represent ? Is it frequency or something else? my question is about the Z- transform. My first question is what the title says. What does 'Z' in Z-transform represent ? Say in Fourier transform, 'w' (omega) represents frequency ? From Fourier transform, I can know what is the strength of the signal at a certain frequency component, right? Then what idea do I get from Z-transform ?

And, in Fourier transform, we can draw a graph with x-axis as the frequency and y-axis as the amplitude of a certain frequency component. Can we do the same for Z-transform, with x-axis as the different 'Z' values and y axis as the Z-transform for a particular Z ?  My textbook says only about equations, but I don't get the significance of 'Z' in  Z-transform. And about the graph thing, I haven't encountered any such graph in any textbook, only the Z-plane.

The Fourier transform doesn't exist for every signal. For example, there is not a Fourier Transform of the signal $ (2)^{n} u(n) $. So, in order to get an useful toll for this signals, you can multiply the signal for a sequence $ r^{-n} $ in order to make it decay fast enough to converge:

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

$$ X(\omega) = \sum_{n=-\infty}^{\infty}{x(n)(re^{j \omega})^{-n}} $$

Then, we can make $ z = re^{j \omega} $. So, in this case, z is a complex value that can be understood as a complex frequency.

It is important to verify each values of $ r $ the sum above converges. These values are called the Region of Convergence (ROC) of the Z transform. The z transform is equal to the Fourier transform when $ r = 1 $, but the Fourier transform will only exist if $ r = 1 $ is inside the ROC.

As $ z $ is a complex variable, you can't plot it in a Cartesian graph as you said. You will need to plot an 3D graph where x-axis and y-axis contain the real and imaginary values of $ z $, and the z-axis the values of the z-transform. I don't know if there is a real utility for a graph like that in any application, and that's why it is not common to see this graph. Normally, the z-plane is enough.


Questions on the generation of R.O.C. curves In my foray into detections, I am trying to understand how best to generate receiver operating characteristic (ROC) curves, and some questions have come up.

In my studies, I have found out that ROC curve generation is simply the plot of true-positive rate (TPR) versus false-positive rate (FPR). Naturally one wants TPR to be 1 and FPR to be 0. I also understand that each point on the ROC curve, (a point with one TPR and FPR coordinate) corresponds to one set of detector 'settings' (threshold(s)). With one setting fixed, we will compare all positive 'detects' to all actual number of true positives, and true negatives, and thus yield TPR and FPR. My questions are as follows:

It seems to me that we must have a labelled data set(s) in order to do this to begin with. That is, we have to manually or otherwise, go through all the data we have, and label them as 'yes' or 'no'. (I realize this might seem obvious but I want to confirm). Is this true?
Somewhat related to the above - if in order to generate a true ROC curve that properly captures the detector performance we have to have 'no' data, then it seems to me that we have to have 'no' data that captures almost every possibility that exists in the universe. Clearly this is not feasible, but then how can one claim with any good faith that a detector truly has a certain false alarm rate? This is what confuses me the most.
Is it unheard of for detectors to have ROC curves that almost never false alarm, but that have normal looking variance on the y-axis? Put another way: Is it unheard of for detectors to have very little variance on the false-positive axis (x-axis) but a lot more variance on the true-positive axis (y-axis), for a fixed SNR?

I think that just about summarizes my questions, thanks.

EDIT Number 1:

I do not think I understand the feedback.

Theoretical aspect:

How can I statistically 'model' every single possibility in order for me to theoretically characterize it to begin with? For example, lets say that I have a detector that detects the presence of a signal of 10 Hz. The input can be anything under the sun. Maybe its 10Hz in noise. Maybe its 11Hz and no noise. Maybe its 10.5Hz on top of a linear chirp, or a bird song. Maybe its brown noise overlayed with human voice, overlayed with a 60Hz power supply. How does one model that, or any of the myriad combinations therein?

I do not know if I can get the theoretical performance of my detector, because quite simply I do not see how thats possible for realistic scenarios where anything can be your input! This goes back to Point#2, whereby, how can one possibly know of every possible input into the detector? I could be wrong on this so please show me how.

Practical aspect:

Therefore in this sense, I suppose my question is: I already have a detector in place. It takes an input vector, and spits out a decision, 'yes' or 'no'. I understand that I have to label my data sets in order to come up with true positive and false positive rates. In this case it seems to me that the only possibility is to come up with a data set that represents just about anything possible under the sun, and test that against your detector.

What am I not getting here?

Question 1) Yes, you must know what the right answer is before-hand in order to judge the receiver's performance.  I am doing this right now on a project.  I generate data, scramble it, modulate it, add noise and carrier offset, etc., then I give it to the receiver.  The bits I started out with are my "yes/no" labels.

Question 2) All such measurements are statistical in nature.  Noise is statistical.  Jitter is statistical.  A priori channel conditions are statistical.  Of course the ROC curve is also statistical.  How could it be otherwise?

Question 3) I don't quite understand this question.  A picture would help.

EDIT: Receivers are not tested in every set of conditions under the sun.  First, the designers rely on everyone obeying the FCC requirements for a given band, so in a cellular band they assume that there won't be any TV stations, or radar chirps, and that all other cellular signals will be at specified frequencies and power levels, with specified maximum carrier offsets.

The designers then model "typical" environments (i.e. channels).  That implies certain transmission losses due to free-space path loss, antenna gains, multi-path, etc.  They then put together various typical and worst-case scenarios and test it under those conditions.  For example, a worst case interference scenario would be having adjacent signals at the maximum allowed power.

If the performance is as expected for those selected conditions, the designers assume, through "interpolation", that it will perform acceptably in the other non-tested environments.  Conditions can always be worse than the worst-case test scenario of course.  That is when the phone drops the call. :-)


What is the Fourier Transform of $\DeclareMathOperator{\rect}{rect} \rect(2Bt)\cos[{\omega}_Ct + k_fm(t_k)t] $? What is the Fourier Transform of $\rect(2Bt)\cos[{\omega}_Ct + k_fm(t_k)t] $?

I got the following as the solution:

$$\DeclareMathOperator{\sinc}{sinc} \frac{1}{2} \frac{1}{2B} \sinc(\frac{\omega+{\omega}_C+k_fm(t_k) +}{4B}) + \frac{1}{2} \frac{1}{2B} \sinc(\frac{\omega-{\omega}_C-k_fm(t_k) +}{4B})$$

However, in the book it is given as:

$$ \frac{1}{2} \sinc(\frac{\omega+{\omega}_C+k_fm(t_k) +}{4B}) + \frac{1}{2}  \sinc(\frac{\omega-{\omega}_C-k_fm(t_k) +}{4B})$$

Wolfram alpha shows this:

https://www.wolframalpha.com/input/?i=Fourier+transform+calculator&assumption=%7B%22F%22%2C+%22FourierTransformCalculator%22%2C+%22transformfunction%22%7D+-%3E%22+rect%282Bt%29cos%28pt%29%22&assumption=%7B%22F%22%2C+%22FourierTransformCalculator%22%2C+%22variable1%22%7D+-%3E%22t%22&assumption=%7B%22F%22%2C+%22FourierTransformCalculator%22%2C+%22variable2%22%7D+-%3E%22w%22

If you define the $\textrm{rect}()$ function according to this definition, then your result is correct.
However, there might be a problem with the definition of the $\textrm{rect}()$ function. If in your book they define that function to have unit area, then the result of the book would be correct.