Does speech volume affect frequencies? I have created a little classifier that listens for certain keywords. Sometimes it will not react to a keyword, but if I shout the word, then the classifier will be able to pick it up.
The classifier uses MFCC as features, which is basically a Short-time Fourier Transform with some extra stuff. This tells me that when I shout, the frequencies that the classifier are detecting are actually different from the non-shout.
So I'm wondering: does the volume of my shout change the frequencies in my speech, or am I somehow enunciating differently when I shout?

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


How to solve a time delay from an FIR filter? I have designed a FIR high-pass filter $\{h_n\}$ with length 257 to filter a signal $\{f_n\}$ with length 512. To convolve signal with filter in time,
I performed an FFT with length 768, padding zeroes at tail of $\{h_n\}$ and $\{f_n\}$, and then I multiplied their frequency responses $\{H_n\}$ and $\{F_n\}$ . To convert it back to time I performed an IFFT producing a filtered signal $\{\tilde{f}_n\}$.
However, $\{\tilde{f}_n\}$ has a time shift comparing to the origin signal, and their amplitudes also have a deviation.
I know there is an 'Overlap-Add Method' to solve this problem, but my calculation resources are limited.
Is there any method to solve such problem using minimal times of FFT/IFFT?

However, {f~n} has a time shift comparing to the origin signal, and their amplitudes also have a deviation.

Any causal system has a non-zero group delay. That's to be expected from systems theory!
So, yes, you've got a shift. And, unless your ${h_n}$ is a linear-phase filter, that delay isn't even constant for different frequencies.
This has nothing inherently to do with your method of convolution – it's just how the math works out.
If your filter is a linear-phase filter, then you can just time-shift the result back by half its length to compensate the group delay.
The usual way, however, is simply to be aware of the math behind all this and to work with the signal knowing that filtering is an operation with a group delay.


How can faux correlated data be generated for testing or training? There are one or two questions here that ask how to assess data correlation. The data tends to be empirical.

Is there some tool or standardised technique for generating correlated data from scratch?  Some way of generating a sequence with an exact amount of known correlation or autocorrelation? It would allow quantitative testing, training and assessment of correlation measurement tools. Or do people just wing it and knock something up Pythonic?

Remember that the power spectral density (PSD) of anything is the Fourier transform of the autocorrelation function of that process.

Hence, if you want a defined autocorrelation function (ACF), do its inverse Fourier transform, and you'll get a PSD that would represent a signal of given ACF.

Now, you can use all the beautiful methods of filter design to come up with a filter that, when applied to a white process (e.g. random numbers from a random number generator) will shape that process's PSD as desired!

So, that's kind of the general look at things: If you just want something like

  "Each sample should be correlated much with the previous one, but only slightly with the one before"

then maybe understanding it as roughly

$$y[n] = 0.9x[n] - 0.1x[n-1]$$

would get you there quicker: $(0.9, 0.1)$ is actually the coefficient vector of a FIR filter that you can apply to your uncorrelated input vector $x$ and be done with it.


Rectangle detection in "real life" images I have seen a lot of posts here and on SO and on various websites, using Matlab, OpenCV, name it, on how to detect rectangles in an image. I have done it myself too, using Matlab. My usual steps were:

Grayscale the image (if not done)
Blur the image (if needed)
Apply an edge detector (e.g. Canny, LoG, etc., tried a lot of them)
Apply the various steps of Hough Transform to find lines
Use some logic in order to tell which line combinations are forming rectangles

I have no problem with these. The problem is that (and with most of other posts I found) it works with relatively simple images, generally synthetics one with few textures, or with high contrast (like a black rectangle surrounding text on a white sheet of paper). Well, maybe that's only me, but I can't make it work with more complex images, where the edgemap might be quite noisy, where I can't afford to manually tune all sort of parameters and thresholds to "make it work".

Here is a dummy project to illustrate the kind of situation I want to address. Imagine that I have a lot of books (dozens) layered and overlapping on a table and I have a camera above, at the right distance to view all of the scene and only the scene. There is no distortion of lens, no blur from the lighting and the resolution is pretty good. The books may be rotated at any angle, they may be of different sizes (though one could assume there is a min and max size), their cover may or may not be different from each others, with possibly images, text, etc. but no colour information is present.

I want to be able to identify the location of every book in that scene. Then, knowing their location, I will be able to extract the appearance of their covers and maybe retrieve their names from a database using these covers (just an example, this part is not important now). Oh, I cannot turn the problem around and use the database to search for the presence of a given book (it would be quite long and I try to constrain this example to rectangle detection).

Now, imagine that, given all the information known about the problem at hand, I conclude that the best solution would be to search for the presence of rectangles in the scene, knowing that books are indeed rectangles and finding all the rectangle could give me a good start. Oh, and I want to do it in real time.

What would you do? I have looked at a some papers from the literature, a lot are doing contour detection to detect horses, mugs, bottles, but they are quite slow to detect the shapes and more I have noticed that often just plain simple rectangles are, well, too simple, not enough distinctive for the algorithms to find.

I am not looking for a ready made answer. I just wonder if anyone has a thought on this. I have not included a list of papers or implementation I tested partly because this question is quite long enough and also I don't want to orient the answers. Remember, this example is made up and what I really want is to detect rectangles. And, disclosure, eventually in my career I assume I will have to detect rectangles (or other shapes) so I was quite ashamed that I wasn't able to do it on real images with the knowledge I have right now   :S

Thanks

I think you would use a 2D Matched Filter. You would convolve your image with a series of rectangles. The peaks in the resulting images would be the location of your books. You could do this quickly in by Fourier transforming your image and using the known function for a rectange in 2D Fourier space (its two sinc functions, multiplied).


Frequency bin of the positive and negative frequency I am using this MATLAB tutorial for Frequency-Domain Linear Regression. There is one part of code they provide, where it is necessary to determine the "frequency bin of the positive and negative frequency", for fft and ifft:

freqbin = 72/12;
freqbins = [freqbin 72-freqbin]+1;
tsfit = zeros(72,1);
tsfit(freqbins) = tsdft(freqbins);
tsfit = ifft(tsfit);
mu = mean(ts);
tsfit = mu+tsfit;

Length of time series is $72$, and $12$ months is considered as one cycle.

How can there be only one frequency bin for positive frequency?
How do we know that it is exactly $72/12 + 1$ ($+1$ is because first bin is for zero frequency)? Is this some formula or what?

In the tutorial it is stated that

  Because the frequencies are spaced at 1/72 and the first bin corresponds to 0 frequency, the correct bin is 72/12+1. This is the frequency bin of the positive frequency. You must also include the frequency bin corresponding to the negative frequency: -1 cycle/12 months. With MATLAB indexing, the frequency bin of the negative frequency is 72-72/12+1.

So, regarding your questions:

It's not that there is only one bin for positive frequencies. The thing is that you are looking for a linear regression model, so you check what frequency is the most significant and then you create what would be the resulting vector if a DFT was made to a discrete sinusoid of that frequency only. That is your model.
The bin number is not any weird formula or anything. You have data for every month, but in the tutorial they are looking for cycles per year. You have 6 years. So $72/12 =6$ and that's where it comes from. (As you said, the $+1$ comes from the fact that there is a bin corresponding to the mean value.)


how to make CIC compensation filter I gotta wrap my head around to design CIC compensation filter

I'm studying by referring these materials:

Altera, "Understanding CIC compensation filters
Hardware Efficient FIR Compensation Filter
For Delta Sigma Modulator Analog to Digital Converters.
Circuits and Systems, 2005. 48th Midwest Symposium on
Saiyu Ren ; Dept. of Electr. Eng., Wright State Univ., Dayton, OH ;    Siferd, R. ;    Blumgold, R. ;    Ewing, R.

first Altera just say "inverse sinc function"
they wouldn't let me know how to realize the function except for using matlab function fir2 and how to get the coefficients of taps

so I found another reference (2) which I mentioned above.

they explain how to decide the coefficients of taps. but I couldn't figure out some equations.

they say

H(z)=a0*z^0+a1*z^-1 + ... + an*z^-n  => H(f)=h(0)+2*sigma(from 1 to (n-1)/2) h(k)*cos(2k*pi*F) where F=f/fs

Here is my first question. How is H(z) converted to H(f). If I substitute z to exp(iwT) it doesn't make sense.

second question I though a0=an=h(1). Is it correct?

finally i used matlab to calculate the coefficients of taps and i got the same coefficients in comparision with a result of the paper, but a simulation gave me a totally different and wrong data

h= [-0.70  2.09 -1.76  0.74 -1.26 0.40  0.34  0.26 0.74
    -0.53  0.24 -0.74  0.04 -0.39 0.05 -0.02  0.28 0.14
     0.41 -0.02  0.36 -0.02  0.41 0.14  0.28 -0.02 0.05 -
     0.38  0.04 -0.74  0.24 -0.53 0.74  0.26  0.34 0.40 -
     1.26  0.74 -1.76  2.09 -0.70] 41 taps

please let me know where I have to modify...

If I violate The copyright issue (because of the content of paper)
please tell me and I will edit or delete this.

Here's my matlab code and plot

clc; clear;
OSR=8
fb=125e6/2

fs=2*fb*OSR
f=[0 14 28 38 45.5 49.5 100 140 170 200 225 250 275  300 325 350 375 400 430 460 490]*10^6;

F=f/fs;
N=41;
L=[1: 1 : (N-1)/2];
S=max(L)

for i=1:S+1
   A(i,1)=1;
end

for i=1:S+1
    for j=2:S+1
    A(i,j)=2*cos(2*pi*(j-1)*F(i));
    end
end

H_sinc=((sin(OSR*pi.*F)./(OSR*sin(pi.*F)))).^2;

H_FIR=1./H_sinc;

H_IFIR(1:6)=H_FIR(1:6);

H_IFIR(1)=1;

H_IFIR(7)=0.1;

H_IFIR(8:21)=0;

H_IFIR=H_IFIR';

A_inv=inv(A);

H=inv(A)*H_IFIR;

H=H';

for i=1:S+1
  H_sol(i)=H(S+2-i);
end

for i=1:20
    H_sol(N-i+1)=H_sol(i);
end
f=[0:100:fs];
F=f/fs;

H_fir1=H_sol(21)+2*H_sol(1)*cos(2*pi.*F) ...
           +2*H_sol(2)*cos(2*pi*2.*F) ...
           +2*H_sol(3)*cos(2*pi*3.*F) ...
           +2*H_sol(4)*cos(2*pi*4.*F) ...
           +2*H_sol(5)*cos(2*pi*5.*F) ...
           +2*H_sol(6)*cos(2*pi*6.*F) ...
           +2*H_sol(7)*cos(2*pi*7.*F) ...
           +2*H_sol(8)*cos(2*pi*8.*F) ...
           +2*H_sol(9)*cos(2*pi*9.*F) ...
           +2*H_sol(10)*cos(2*pi*10.*F) ...
           +2*H_sol(11)*cos(2*pi*11.*F) ...
           +2*H_sol(12)*cos(2*pi*12.*F) ...
           +2*H_sol(13)*cos(2*pi*13.*F) ...
           +2*H_sol(14)*cos(2*pi*14.*F) ...
           +2*H_sol(15)*cos(2*pi*15.*F) ...
           +2*H_sol(16)*cos(2*pi*16.*F) ...
           +2*H_sol(17)*cos(2*pi*17.*F) ...
           +2*H_sol(18)*cos(2*pi*18.*F) ...
           +2*H_sol(19)*cos(2*pi*19.*F) ...
           +2*H_sol(20)*cos(2*pi*20.*F);

figure(1), semilogx(f,db(H_fir1))
% hold on
grid on

please help me to solve this problem...

I have a simple compensation approach that I've used to implement a reasonable estimate of an inverse $\textrm{sinc}$ for use as a CIC compensator as well as other inverse $\textrm{sinc}$ applications. This approach makes use that the $\textrm{sinc}$ function in the passband can be reasonably approximated by a weighted cosine function. Therefore a raised cosine function in frequency, meaning $b-(\alpha)cos(\omega)$ can be used in cascade with the CIC filter as compensation. This function is adjusted by changing the cosine weight $\alpha$ while keeping $b=1+\alpha$, such that the mean squared error over the passband of interest of the cascaded result is minimized.

Once the minimized passband is established, (you can do this with least squared minimization techniques, but I have done it quickly by evaluating passband error for a least squared miminum while doing a binary search on $\alpha$ values, which quickly converges to an $\alpha$ value that minimizes the error for a given passband), the coefficients of the filter, which is the filter impulse response, are as follows:

Coeff 1, 3: $-\alpha/2$

Coeff 2: $b= 1 + \alpha$

It is that easy!  The filter is a simple 3 tap FIR with coefficients [$(-\alpha/2)$   $(1+\alpha)$   $(-\alpha/2)$]. Since the filter is symmetric, it can also be implemented with just two multipliers as shown in the figure. Symmetric filters also have the nice benefit that they are linear phase.

See the figure below for an example of using this CIC compensator in an x8 CIC Intepolator.  In use, the signal is passed through the 3 tap compensator at the lower rate, and then fed into the x8 CIC Interpolator with an output at 8x the sample rate.  The cascaded response shows the significant improvement that can be achieved with this simple compensator.

What should now be clearer from this figure is that this compensation approach will slightly decrease overall rejection in portions of the stop band. (For this specific example the decrease was on the order of 1.6 dB). However, since the unity gain positions of the compensator will be aligned with the null locations in the CIC response, this reduction in rejection is not in those rejection regions that matter most for CIC rate conversion.

Here is a zoom in of the passband showing the excellent flatness that can be achieved with three taps!

Note on solving for the FIR coefficients from the first figure above.

The coefficients being the impulse response of the filter are determined from the inverse DFT of the frequency response. They can be quickly derived from the generic equation for the frequency response for an FIR filter in terms of the FIR coefficients as:

$H(\omega)= \sum_{n=0}^{N-1}c_ne^{-jn\omega}$

where $c_n$ are each of the $N$ coefficients, and $\omega$ is the digital frequency domain from $0$ to $2\pi$ corresponding to 0 Hz to the sampling rate.

As for the case of our 3 tap symmetric FIR this is

$H(\omega)= -\frac{\alpha}{2}+ be^{-j\omega}-\frac{\alpha}{2}e^{-j2\omega}$

$ = e^{-j\omega}(-\frac{\alpha}{2}e^{+j\omega}+b-\frac{\alpha}{2}e^{-j\omega})$

$=e^{-j\omega}(b-\alpha(\frac{e^{+j\omega}+e^{-j\omega}}{2}))$

$=e^{-j\omega}(b-\alpha cos(\omega))$

which is exactly the desired frequency response as given in the first plot (the first term $e^{-j\omega}$ is just the required delay for the filter to be causal but does not affect the magnitude response (and shows how the phase is linear versus $\omega$ ...the same reason all symmetric FIR filters are linear phase- the response for any symmetric FIR can all be described in terms of cosines with a linear phase delay pulled out of the equation!).

It may be more intuitive to some imagine this operation when going from the time domain to the frequency domain if you are familiar with the FT of a cosine wave: A raised cosine wave in the time domain (meaning a cosine with a DC offset) has a Fourier tranform with 3 impulse components in the frequency domain; a DC term, and a positive and negative frequency term. Likewise, the same waveform in the frequency domain will have three impulse components in the time domain, and these components are the impulse reponse by definition, which are therefore the filter coefficients.


How to understand the de-correlation property of DCT? What does de-correlation mean? Here are three problems that I can't understand:

Usually what does de-correlation mean in image-processing?
What's the benefit of de-correlation when it comes to image splicing?
Is there some relationship among neighbouring coefficients? And I read a paper which says coefficients difference matrix(got by coefficients subtracting its neighbor) can eliminate the relationship of the coefficients，How to understand this?

I would really appreciate if someone can help me!

The term "decorrelation" generally accounts for processing that reduce autocorrelation within single, or cross-correlation within a group of  signals or  images. Doing that, it should preserve important features in the  data.

In other terms, images often contains simpler "objects" (in terms of  morphology: bumps, edges, textures) that are mixed, or distorted. Generally, unmixing them via decorrelation somehow simplify further processing.

Generally yes, unless you have uncorrelated noise in your images. A classical and simple model is an order-one Markov or autoregressive process. Basically a pixel somehow  depends (with a parameter $\rho$) on the past pixel, with uncertainty. If $\rho$ is close to $0$, you are already decorrelated. If $\rho$ is close to $1$ (typically $0.90-0.95$), the  data is quite correlated.

The DCT was meant for diagonalizing the  resulting autocovariance matrices with  Toeplitz structure, to give fast estimates of their eigenvectors. Transforming a fat autocovariance matrix into a thin, close to a diagonal, matrix is an instance of decorrelation.


Summations in Z-Transform I'm currently working on a problem that involves a Z-Transform. Basically, the essence of the problem is that if:
\begin{equation}
H\left(z\right) = \sum_{n=0}^{N-1}h\left(n\right)z^{-n}
\end{equation}
find:
\begin{equation}
Y\left(z\right)\:=\:H\left(z\right)\cdot H^*\left(\frac{1}{z^*}\right)
\end{equation}
The asterisk represents a complex conjugate.
So, I know that
\begin{equation}
H^*\left(\frac{1}{z^*}\right)\:=\:\sum _{n=0}^{N-1}h^*\left(n\right)z^n
\end{equation}
Am stuck at this part now:
\begin{equation}
Y\left(z\right)=\left(\sum_{n=0}^{N-1}h\left(n\right)z^{-n}\right)\left(\sum_{n=0}^{N-1}h^*\left(n\right)z^n\right)
\end{equation}
I'm stuck at this part and don't know what to do from here. Thank you for your help in advance.

This may help: the result is the cascade of an FIR filter with coefficients $h(n)$ with the reverse conjugate filter with coefficients $h^*(N-n)$. If we can limit $H(z)$ to be the unit circle (for computing the frequency response), then
$$H(z^*)|_{z=e^{j\omega}} = H(z^{-1})|_{z=e^{j\omega}}$$
since
for $z = e^{j\omega}$, $z^* = e^{-j\omega}$
And therefore
$$H^*(1/z)|_{z=e^{j\omega}}= H^*(z^*)|_{z=e^{j\omega}}$$
And
$$H^*(1/z^*)|_{z=e^{j\omega}}= H^*(z)|_{z=e^{j\omega}}$$
So
$$Y(z)|_{z=e^{j\omega}} = H(z)H^*(z)|_{z=e^{j\omega}} = |H(z)|^2|_{z=e^{j\omega}}$$

Also from the link the OP provided to me in the comments below, it is clear that the problem is dealing with the two sided Z-transform and the full problem has a solution given as
$$H(z)H^*(1/z^*)\Phi_{xx}(z) = \frac{\sigma_v^2}{(1-az^{-1})(1-a^*z)}$$
And we see from the earlier solution in the link that
$$H(z) = \frac{1}{1-az^{-1}}$$
Therefore
$$H^*(z) = \frac{1}{1-a^*(1/z^*)}$$
And
$$H^*(1/z) = \frac{1}{1-a^*(z^*)}$$
And
$$H^*(1/z^*) = \frac{1}{1-a^*z}$$
And
$$H(z)H^*(1/z^*)=\frac{1}{(1-az^{-1})(1-a^*z)}$$
Note confirming my earlier result, when we restrict z to be the unit circle (for determining the magnitude response):
$$H^*(z)|_{z=e^{j\omega}} = \frac{1}{1-a^*(1/z^*)}|_{z=e^{j\omega}}= \frac{1}{1-a^*(1/e^{-j\omega})} = \frac{1}{1-a^*(e^{j\omega})} = \frac{1}{1-a^*z}|_{z=e^{j\omega}}$$
Also to avoid any confusion, earlier in the linked solutions where it says "from 7... we arrive at (2.55) as follows" it gives a result that indicates
$$\Phi_{xx}(z) = \Phi^*_{xx}(1/z^*)$$
But this is not generally true for any z transform, (we cannot conclude that $H(z) = H^*(1/z^*)$ for all z) but it is in this case specific to the autocorrelation function $\phi_{xx}(k)$ since that function is Hermitian:
$$\phi^*_{xx}(k) = \phi_{xx}(-k)$$


Tie breaking in Otsu method I had a doubt regarding the selection of minimum mixed variance, if a situation arises that same minimum mixed variance occurs two times. Is there any way to break this tie?

You need to add another criteria for the tie breaking.

Maybe how coherent each threshold area is after the we apply the threshold.

So for example lets say you have a tie between thresholds T1 and T2. You can apply T1 and than find the number of connected component C1 it gave you than apply T2 and find the number of connected component C2 it gave you and than decide to take the threshold that gave less connected components.


Using SURF algorithm to match objects on MATLAB I am trying to find the difference between two images, using Matlab. The classic built in function that Matlab provides for this is because the two images don't have the same dimensions (The objects in the images are the same, but in the second image other objects are introduced).
i need to use a metric or function can calculate the number of matching feature between two images and give me a value to decide according to threshold ..

source = imread('source.png');
target = imread('target.png');
source = rgb2gray(source);
target = rgb2gray(target);
  sourcePoints=detectSURFFeatures(source,'MetricThreshold',100.0,'NumOctaves',1,'N   umScaleLevels',6);
 targetPoints=detectSURFFeatures(target,'MetricThreshold',100.0,'NumOctaves',1,'NumScaleLevels',6);
     [sourceFeatures,sourcePoints]=extractFeatures(source,sourcePoints,'SURFSize',64);
     [targetFeatures,targetPoints]=extractFeatures(target,targetPoints,'SURFSize',64);

  boxPairs = matchFeatures(sourceFeatures, targetFeatures);

  matchedSourcePoints = sourcePoints(boxPairs(:, 1), :);
  matchedTargetPoints = targetPoints(boxPairs(:, 2), :);

  % IS there any metric can used to make accurate decision if there is a match or not instead of percentage because its not accurate as required
  numPairs = length(boxpairs); %the number of pairs
  percentage  = numPairs/100;

     if percentage >= 0.40
     disp('We have this');
       else
     disp('We do not have this');
     disp(percentage);
    end

IS there any metric can used to make accurate decision if there is a match or not instead of percentage because its not accurate as required

People generally use Euclidean distances, just like you do and it works as it is, pretty well. However, if you require further robustness, I propose the following:

Compute the matching both ways: image1->image2 and image2->image and then retrieve only the agreeing matches. This introduces some consensus.
Assuming you have planar patches, you run a homography using RANSAC. The good thing here is that RANSAC is robust to outliers, given a predifined tolerance. This way prune away a lot of matches, and you will be left with a meaningful set of inliers.

Using these two together should really give you robust and repeatable matches.


Is there a process for deriving special cases of convolution? For my signals and systems homework, I was met with the following problem:

Consider a CT system with impulse response given by

$$h(t) = e^{-|t|}$$

and a step input $x(t) = u(t)$. Using the deﬁnition of convolution, determine an expression for the output $y(t) = x(t)∗h(t)$.

I approached this problem using the standard method taught in my class. I set up the following integral:

$$\int_{-\infty}^{+\infty}h(\tau)x(t-\tau)d\tau$$

While evaluating the integral, I got stuck at the end as the integral approached infinity. I decided to rearrange $h(t)$ by expressing it in a different form. The form I then got it in was:

$$h(t) = e^{-t}u(t) + e^tu(-t)$$

Which then allowed me to match it to the following convolution on a table given with my homework.

From doing this problem it seems to me that there are special cases for convolutions. Ones that which the standard method of evaluation doesn't work.

Is there a process to derive these special cases? If so can an example be given? And why wouldn't the standard process work for this example?

Well, indeed there are special cases for convolutions but yours is quite straightforward. You have to consider two cases: $t<0$ and $t>0$.

If you sketch the convolution process for these two cases, by sliding your time shifted and time reversed input on the same axis as your impulse response, you will end up with the following:

Can you handle the integrals shown?


Why does diagonal loading of a covariance matrix make an adaptive beamformer more robust in the case of a perturbed array? It has been shown that 'diagonal loading' a covariance matrix derived for an adaptive beamformer can improve robustness of the beamformer when the antenna array is perturbed, albeit at the expense of background noise power.

The question is why does this work? Is the diagonal loading stopping the adaptive filter from reaching a local minima? Or is there some other mechanism at work here?

I had worked on an array processing problem where I had used diagonal loading of the measurement covariance matrix. But I had used diagonal loading as a solution to what I thought was a numerical issue with my eigen values being too small. Since I had to invert the covariance matrix with small eigen values, I had numerical issues with the result.
Diagonal loading of the covariance matrix helped me get rid of the issue as the eigen values were now bounded by the diagonal loading value.

I am not sure, but one place where diagonal loading helps is in numerical stability of inverting the covariance matrix which is most likely encountered in adaptive beamforming.


What is the definition of half duplex?Is this system,$A$ give $x$ to $B$,then $B$ give $x$ to $C$,called half duplex? I saw the definition of half duplex in wiki,and it said

half-duplex (HDX) system provides communication in both directions, but only one direction at a time (not simultaneously). Typically, once a party begins receiving a signal, it must wait for the transmitter to stop transmitting, before replying.

However,i saw a paper: http://jultika.oulu.fi/files/nbnfi-fe2019040110636.pdf ,its system model is a transmitter transmit the $x_1$ and $x_2$ information to the user1(U1),then U1 transmit $x_2$ to the user2(U2),and it said this is half-duplex.

So in fact,the definition of half duplex is that device can transmit and receive,but no need to transmit and receive at the same time?

I mean,$A$ give $x$ to $B$,then $B$ give $x$ to $C$,this is also called half-duplex?I thought this is called simplex duplex before,because i thought the half-duplex is $A$ give $x$ to $B$,then $B$ give $x$ back to $A$.

So i want to ask that $A$ give $x$ to $B$,then $B$ give $x$ to $C$,this is also half-duplex?

The receiver/transmitter pair in B's possession is operating in half-duplex mode.

As you say "$A$ give $x$ to $B$, then $B$ give $x$ to $C$" which specifies that while $B$ ls listening to what $A$ is telling it (the $x$ that is being given to $B$), the transmitter in $B$'s possession must remain silent until $B$ has completely received $x$. It is only then that $B$ can start up his transmitter and send $x$ to $C$. So, $B$ cannot receive and transmit at the same time, that is, $B$ is operating in a half-duplex mode. If $A$ were to send $y$ to $B$ while $B$ is sending $x$ to $C$, $B$ just ignore the transmission from $A$: $B$'s transmitter is active and so its receiver must be in-operational till $B$ is done talking to $C$.


Why real part of FFT converts image into rotation + original? I have read this image:

taken its FFT (2D) and then Inverse FFT to get exactly the image back. Code is provided for reference:

imfft = fft2(photographer);
im = uint8(ifft2(imfft));

imshow(im); %Output is same image

But when I change the fourier and take only the real part,

imfft = real(fft2(photographer));
im = uint8(ifft2(imfft));
imshow(im);

I get an image like this (note that size change is irrelevant and only due to saving it from Matlab figure handler):

Can someone explain me the theory (maths) behind it? Thanks

Let's say your image is given by $I(x,y)$. Then its Fourier transform is given by
$$
I^f(\omega_x,\omega_y) = \int_x\int_yI(x,y)e^{j\omega_xx}e^{j\omega_yy}dxdy
$$

Now you take the real part and perform the inverse:

\begin{align}
I_m(\alpha,\beta) &= \int_{\omega_x}\int_{\omega_y}\Re\left\{I^f(\omega_x,\omega_y)\right\} e^{j\omega_x\alpha}e^{j\omega_y\beta}d\omega_xd\omega_y
\\
&= \int_{\omega_x}\int_{\omega_y}\Re\left\{\int_x\int_yI(x,y)e^{j\omega_xx}e^{j\omega_yy}dxdy\right\} e^{j\omega_x\alpha}e^{j\omega_y\beta}d\omega_xd\omega_y\\
&= \int_x\int_yI(x,y)\int_{\omega_x}\int_{\omega_y}\Re\left\{e^{j\omega_xx}e^{j\omega_yy}\right\} e^{j\omega_x\alpha}e^{j\omega_y\beta}d\omega_xd\omega_ydxdy
\end{align}

You can clearly see that the inner integral is the 2D Fourier transform of
$$
\cos(\omega_xx)\cos(\omega_yy) + \sin(\omega_xx)\sin(\omega_yy)
$$
which is
$$
\frac{1}{2}\left[\delta(x-\alpha)\delta(y-\beta) + \delta(x+\alpha)\delta(y+\beta)\right]
$$

Substituting the result to $I_m$ yields
$$
I_m(x,y) = \frac{1}{2}\left[I(x,y)+I(-x,-y)\right]
$$

Of course in your case $x,y>0$, however the discrete Fourier transform assumes your signal is $N$-periodic and you get
$$
I_m(x,y) = \frac{1}{2}\left[I(x,y)+I(N-x,M-y)\right]
$$
where $N,M$ are the dimensions of your image. I think you can see now why you got that result.


Generate binary signal and convert to integer in Matlab I am going to test an algorithm and need a test signals.
My first test signal is a binary signal of 1 or -1. I was thinking about how to implement such a sequence in matlab and wrote this:
binary_data = rand(1,M) > 0.5;
b = 2.*binary_data - 1;

My next step is:
dataInMatrix = reshape(b,length(b)/k,k); % reshape data into binary k-arrays, k=log2(M), M -size of signal constelation
dataSymbolsIn = bi2de(dataInMatrix);

error:
Expected input number 1, binary vector b, to be nonnegative.
How to convert such binary sequence to integer if i need create symbols?
EDIT 1.
Guys, I have generated a test signal and it look like this:

i generate a sequence 1 or -1
use a raised cosine filter

How can I get a blue lines...show samples ( digital signal)?

@Ali23. If by "How can I get a blue lines?" you mean how to generate samples of the your red curve, try the following:
m = [1, -1, -1, 1 1, -1]; m = [m m m m];
N = length(m);
M = 10; mup=zeros(1,NM); mup(1:M:NM)=m;  % oversample by M
h = hamming(M);                     % blip pulse of width M
x = filter(h,1,mup);     % convolve pulse shape with 'm' data
t = 1/M:1/M:length(x)/M;
figure(1), subplot(2,1,1), plot(0:0.1:0.9,h)
subplot(2,1,2), plot(t,x)
In the above I used a simple hamming sequence rather than the impulse response of a RRC filter. (You should use the impulse response of a RRC filter.) The key point here is that you have to upsample (zero stuffing) the 'm' sequence before you apply it to the filter.


Is this system an FIR filter or an IIR filter? Consider this system:

 y[n]=y[n-1] + x[n]

It can be verified easily that this system has a finite impulse response (by putting x[n]= delta[n]  i.e impulse, and evaluating impulse response h[n]=y[n]). Yet for a unit step input U[n], this system is an unstable system.  On wikipedia it is mentioned that FIR filters are always stable

Now my question is , is this an FIR filter or an IIR filter?

It is an IIR filter.  I think you misunderstand what is meant by "finite impulse response".  It does not mean finite amplitude, it means finite time- i.e. that the impulse response eventually becomes zeros.

The impulse response of the filter is a step function, which is infinite in time.

x[-1] = 0
x[ 0] = 1
x[ 1] = 0 ...

y[-1] = 0
y[ 0] = 1
y[ 1] = 1
y[ 2] = 1 ...


Showing error energy goes to zero Let
$$\hat{x}[k] = \frac{1}{2\pi}\int_{-W}^{W}X(e^{j\omega})e^{j\omega k}d\omega,\label{ift}\tag1$$
where
$$X(e^{j\omega}) = \sum_{n=-\infty}^{+\infty} x[n]e^{-j\omega n}\label{dft}\tag2$$
Also,
$$d[k] = x[k] - \hat{x}[k]\text.\label{error}\tag3$$

Show that energy of $d[k]$ goes to zero as $W\to \pi$. That is,
$$\lim_{W \to \pi} \sum_{k = -\infty}^{+\infty}|d[k]|^2 = 0\text.$$

My Attempt
Main problem for me is removing magnitude in some way. Writing $|d[k]|^2 = d[k]d^{*}[k]$ complicates equations. Also if we use Parseval's identity, the magnitude is still present. I don't know if we are allowed to move limit inside the summation or not. Even assuming it's possible, I couldn't found a way for calculating $\lim_{W \to \pi} |d[k]|^2$.

HINT: The error energy can be written as
$$\begin{align}\sum_k\big|d[k]\big|^2&=\sum_k\left|x[k]-\frac{1}{2\pi}\int_{-W}^W\sum_nx[n]e^{-jn\omega}e^{jk\omega}d\omega\right|^2\\&=\sum_k\left|x[k]-\sum_nx[n]\frac{1}{2\pi}\int_{-W}^We^{j(k-n)\omega}d\omega\right|^2\tag{1}\end{align}$$
Now compute the integral
$$I(k-n,W)=\frac{1}{2\pi}\int_{-W}^We^{j(k-n)\omega}d\omega\tag{2}$$
and show that $\lim_{W\to\pi}I(k-n,W)=\delta[k-n]$. This reduces the second sum in $(1)$ to a single element and the result follows.


What are the main advantages and disadvantages of the lifting scheme? Can anybody tell me the advantages of the lifting scheme and why it is better then a normal filter bank?

Unlike the convolution, the lifting scheme can

map integers to integers (CDF 5/3 transform in case of JPEG 2000),
asymptotically reduces the computational complexity by a factor two,
can be computed in-place,
allows to simply treat signal boundaries (periodic symmetric extension in case of JPEG 2000).


Minimum Channel Bandwidth for PCM Please check if my answer to GATE IN 2010 question 11.20 is correct/can be improved upon.

So, every $\frac{1}{f_s}$ sec, a sample is taken. The signal amplitude is quantized into $2^n$ levels such that each sample can be encoded using $n$ bits. The question (I think) asks what is the minimum bandwidth a channel needs in order to transmit these bits.

We have $\frac{1}{f_s}$ sec time and $n$ bits to send. So the bit rate would be $nf_s$ bits per second. Since each "bit" is sent as a pulse across the channel, the pulse rate $f_p = nf_s$. So by the Shannon–Hartley theorem we can conclude,

$$nf_s \leq 2B$$
$$B \geq \frac{nf_s}{2}$$

So the answer is (c)

The question asks what the bandwidth needed for "faithful reconstruction" is, I guess they mean that if you don't have at least this large a bandwidth errors will creep in and you can't reconstruct the signal from the n bit code at the receiver end?

I don't think the exam question is stated properly, unless there is additional context that the student is supposed to know.

First, the bandwidth depends on the encoding. The question seems to assume binary encoding, but PCM does not specify this.

Second, the bandwidth depends on the pulse shape. Usually one assumes rectangular unless otherwise specified, but this information should bave been included in the question.

In general, the PCM signal can be written as a pulse train:

$$s_{\text{PCM}}(t) = \sum_{k=-\infty}^\infty a_k p(t-kT_p),$$

where $T_p$ is the pulse rate and $a_k$ is the encoded PCM data. In general, the bandwidth of $s_{\text{PCM}}(t)$ is equal to the bandwidth of $p(t)$.

As an example, let's say we want to transmit the bits $1010$ using a bipolar, non-return to zero binary encoding with $a_k \in \lbrace 1, -1 \rbrace$. The pulse shape is rectangular with duration $T_p = 1/(nf_s)$. Assuming $n$ starts at $n=1$, the PCM signal can be written as:

$$s_{\text{PCM}}(t) = p(t-T_p)-p(t-2T_p)+p(t-3T_p)-p(t-4T_p).$$

The problem is, what is the bandwidth of this signal? Theoretically, it's equal to the bandwidth of the rectangular pulse $p(t)$ with duration $T_p$. However, there are several ways to measure that bandwidth.

For instance, if you use the "absolute bandwidth" definition, then the bandwidth is infinite. There are ten or more engineering definitions of bandwidth, and each gives you a different number. Note, however, that transmitting $s_{\text{PCM}}(t)$ over any channel with less than infinite bandwidth will result in some amount of distortion.

Now let's say that $p(t)$ is a sinc pulse of appropriate width. The expression for the signal doesn't change. Theoretically, now the signal has infinite time duration, but its bandwidth is finite and equal to $1/(2T_p)=nf_s/2$.


Unexpected peaks in spectral analysis I am doing the spectral analysis of some antenna recorded signals by computing PSD (logarithmic absolute value of the FFT) and the spectrogram of the signal. The problem is that I get some unexpected spikes in the signal that appear through all the signal, as it can be seen in the next figure:

The two peaks that appear are the spikes I am speaking about, and they are the same in all the signal. The center spike is not located in $0\textrm{ Hz}$, so it is not a DC offset of the signal. This is giving me some problems, so I would like to know what can they be or how can they be removed (maybe I am not computing the spectral analysis correctly).

I have also performed Welch PSD estimations, and the peaks still appear. I have been told that it can be some internal interference form the electronics recording the data or some issue with ADCs.

The large spike at ~9MHz could be some form of interference from electronics near the receiver or internal to the receiver (maybe from using a poorly conditioned power source).  My best guess for the spike at DC is a DC bias from quantization.  If the digitizer truncates instead rounds, the resulting digital signal will have a non-zero mean and therefore a DC offset.  We had this happen in an FPGA downcoverter implementation at a previous job.  Once they changed the FPGA to round, the DC offset went away.  If you have ever used an Ettus USRP radio, you see the same DC offset in many of those boxes.


Two question about RBJ Peak filter I'm new into DSP and Filters in general. Its the first step ahead I do from learning and using a 1-pole IIR basic Low Pass Filter.

I'm about to implement for my VST audio application (using C++ and IPlug Framework) a RBJ's Peak filter.

By steps (from that document), the formula that I'll use within my Process() function will be y[n], thus (as code):

inline double Process(double input) {
    mBuffer += (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2] - (a1/a0)*y[n-1] - (a2/a0)*y[n-2];
    return mBuffer;
}

Since I'll use peakingEQ, I need to calculate these coefficients:

b0 =   1 + alpha*A
b1 =  -2*cos(w0)
b2 =   1 - alpha*A
a0 =   1 + alpha/A
a1 =  -2*cos(w0)
a2 =   1 - alpha/A

which all depends by user defined parameters and intermediate variables.

My questions are:

What's the difference between using Q (thus sin(w0)/(2*Q)) or BW  (thus sin(w0)*sinh( ln(2)/2 * BW * w0/sin(w0))) on calculating alpha? Isn't Bandwidth = Q when reasoing as Peak filter?
What's the meaning of H(s) function (which is (s^2 + s*(A/Q) + 1) / (s^2 + s/(A*Q) + 1) for the peaking filter)? I don't see where to "insert" it in the formulas.

I hope you can help me to undestand these two dubts.

Maybe this thread helps you in this Q - Prewarping both resonant frequency $f_0$ and bandwidth (or $Q$) when using bilinear transform
H(s) is the s-plane transfer function

Analog prototype of RBJ peak filter (Octave):

K  = 10^(gain/40);
w0 = 2*pi*fc;
% H(s) = (s^2 + (K w0/Q)*s + w0^2) / (s^2 + w0 s/(K*Q) + w0^2)
BC = [1 K*w0/Q w0^2];
AC = [1 w0/(K*Q) w0^2];


Output of lowpass filter with damped sine wave input The random process $$Y(t)=\cos(\omega_0t)\cos(\omega_0t+\pi N(t))$$where $N(t)$ is a Poisson process of parameter $\lambda$ enters a lowpass filter with transfer function $$H(j\omega) =
\left\{
 \begin{array}{ll}
  1  & \mbox{if } |\omega| \leq \frac{\omega_0}{2} \\
  0  & \mbox{if } |\omega| \geq \frac{\omega_0}{2}
 \end{array}
\right.$$ If $Z(t)$ is the output of the filter, find the expectation and variance of both $Y(t)$ and $Z(t)$.

So, first we have that:

$$\mathbb{E}[Y(t)]=\cos(\omega_0t)\mathbb{E}[\cos(\omega_0t+\pi N(t))]=\cos(\omega_0t)\mathbb{E}[\cos(\omega_0t)\cos(\pi N(t))-\sin(\omega_0t)\sin(\pi N(t))]]=\cos^2(\omega_0t)\mathbb{E}[\cos(\pi N(t))]$$

Doing that expectation by definition (with the summation) we get that $$\mathbb{E}[Y(t)]=\cos^2(\omega_0t)e^{-2\lambda t}$$

With a similar procedure we find that

$$\mathbb{V}[Y(t)]=\cos^4(\omega_0t)\cdot(1-e^{-4\lambda t})$$

First I tried to find $\mathbb{E}[Z(t)]=h(t)*\mathbb{E}[Y(t)]$

So using the property that states that the transform of $x(t)y(t)$ is $\frac{1}{2\pi}X(j\omega)Y(j\omega)$ we get that

$$\mathcal{F}(\cos^2(\omega_0t))=\frac{\pi}{2}\delta(\omega+2\omega_0)+\frac{\pi}{2}\delta(\omega-2\omega_0)+\pi\delta(\omega)$$

But it is known that the transform of the exponential multiplied by the step function (the random process involved is valid only for $t>0$ by the definition of the Poisson process probability function) is $\frac{1}{2\lambda+j\omega}$ so

$$\mathcal{F}(y(t))=Y(j\omega)=\frac14\frac{1}{2\lambda+j(\omega+2\omega_0)}+\frac14\frac{1}{2\lambda+j(\omega-2\omega_0)}+\frac12\frac{1}{2\lambda+j\omega}$$

So what I need would be to antitransform $H(j\omega)Y(j\omega)$, and now I'm stuck. How could this be done? And regarding $\mathbb{V}[Z(t)]$, I have no idea of how to even begin to solve it.

I've finally figured this out thanks to a friend of mine. The thing is that $N(t)$ can be even or odd, and in both cases, the signal we get is easy to deal with.

We can calculate the desired expectation this way:

$$\mathbb{E}[Z(t)]=\mathbb{E}[h(t)*Y(t)]=\mathbb{E}[h(t)*Y(t)|N(t) =2n]\mathbb{P}[N(t)=2n]+\mathbb{E}[h(t)*Y(t)|N(t) =2n+1]\mathbb{P}[N(t)=2n+1]$$

But when $N(t)$ is even or odd, $Y(t)$ becomes simple:

$$Y(t) =
\left\{
 \begin{array}{ll}
  \cos^2(\omega_0t)  & \mbox{if } N(t)\mbox{ is even } \\
  -\cos(\omega_0t)\sin(\omega_0t)    & \mbox{if } N(t)\mbox{ is odd }
 \end{array}
\right.$$

So

$$\mathbb{E}[Z(t)]=\mathbb{P}[N(t)=2n]\cdot h(t)*\cos^2(\omega_0t)-\mathbb{P}[N(t)=2n+1]\cdot h(t)*(\cos(\omega_0t)\sin(\omega_0t))$$

And those two probabilities are easy to calculate, as well as the two convolutions left.


Relationship between real and imaginary part of a real-valued and causal system I have one question about the real part of a real-valued and causal system with the imaginary part of its Fourier transform given by

$$\textrm{Im}\big\{X(e^{j\omega})\big\}=3\sin(2\omega)-2\sin(3\omega)$$

How can I obtain the real part of the Fourier transform?

i have simulated this by octave here:

>> w=[1:100];
>> Im_X=sin(2.*w)+sin(3.*w);
>> plot (Im_X);

Is it possible to solve in by the Octave/Matlab?

I suppose you've learned about the Hilbert transform relationship between the real part and the imaginary part of the Fourier transform of a real-valued and causal signal. So you can obtain the real part from the imaginary part using this relationship. Since the imaginary part is just composed of sinusoids, the real part can be obtained very easily (note that the Hilbert transform of a sine is a negative cosine).


Choosing samples per symbol for modem pulse One of the parameters in my DSP library for pulse design is "samples per symbol", and I would like to get some advice about choosing this parameter when designing a modem.

The smaller the pulse width, the more pulses you can transmit per second, and so the higher the data rate. Therefore, to maximise data rate I should minimise the pulse width. But just how small can I make the pulse?

The closest I came to finding an answer to this is this, but does not directly address the choice of pulse width in samples.

The pulse rate is ultimately limited by the channel bandwidth. Assuming you're using raised cosine pulses with rolloff factor $\beta$, then $$B = \frac{(1+\beta)R_p}{2},$$ where $R_p$ is the pulse rate.

The number of samples per pulse is, in all libraries I know, independent of the pulse rate. Strictly speaking, you don't need more than $2B$ samples per second. However, some algorithms work better with more samples per pulse, and any signal plots you create will be clearer. In addition, normally you don't want your DAC to be working close to the theoretical limit.

In my experience (mainly with GNU Radio and USRP radios), 4 samples per pulse works very well with root raised-cosine pulses.

To further clarify, let us work through an example. Assume we want to transmit at a pulse rate $R_p=10,000$ pulses per second, with $\beta=0.5$. We know that we'll need a bandwidth of $7,500$ Hz. If we set our system to use 4 samples per pulse, then we'll have a sampling rate of $40,000$ samples per second. This is the amount of samples that our DSP system will need to process per second. If, instead, we go with 6 samples per pulse, then the sampling rate will increase to $60,000$ samples per second and the computational load of the system will increase accordingly.


Inverse Fourier transform of complex exponential with frequency dependent shift In the case of a constant delay $\tau$, we have the following equality:
$$\begin{align}\mathcal{F^{-1}}\left\{e^{-j\omega \tau}\right\}=\delta(t-\tau)\end{align}$$
If the delay is frequency dependent $\tau(\omega)$, can ${F^{-1}}\left\{e^{-j\omega \tau(f)}\right\}$ be expressed as a sum of diracs?

In the general case: no.
For the (inverse) Fourier transform of a function to be composed of a sum of countable diracs (i.e. to be discrete), the function needs to be periodic.
Your $e^{-j\omega\tau(f)}$ is not periodic (it can be periodic, if $\omega\tau(f)$ happen to be periodic with period rationally related to $2\pi$, but that's a pretty special case). Therefore, its inverse Fourier transform is continuous, and can't be represented by a sum of diracs.


normalized Laplacian of Gaussian Laplacian of Gaussian formula for 2d case is
$$\operatorname{LoG}(x,y) = \frac{1}{\pi\sigma^4}\left(\frac{x^2+y^2}{2\sigma^2} - 1\right)e^{-\frac{x^2+y^2}{2\sigma^2}},$$ in scale-space related processing of digital images, to make the Laplacian of Gaussian operator invariant to scales, it is always said to normalize $LoG$ by multiplying $\sigma^2$, that is
$$\operatorname{LoG}_\text{normalized}(x,y) = \sigma^2\cdot \operatorname{LoG}(x,y) = \frac{1}{\pi\sigma^2}\left(\frac{x^2+y^2}{2\sigma^2} - 1\right)e^{-\frac{x^2+y^2}{2\sigma^2}}.$$ I wonder why multiply by $\sigma^2$ not $\sigma^4$ or anything else?

Laplacian response decays as $\sigma$ increases, and it is the second Gaussian derivative so it is multiplied by $\sigma^2$.

LOG is defined as $\bigtriangledown^2G$ so the scale normalized LOG would be $\sigma^2\bigtriangledown^2G$. You need to get rid of the scaling factor of the Gaussian which is $\sigma^2$.


Baseband representation of a LFM chirp signal I need to generate a chirp signal between $f_{start} = -f_{limit}$ and $f_{stop} = f_{limit}$, to compare it to a FrFT generated signal.

The FrFT signal of order $\alpha = 0$ is PAM modulated data signal, so is in baseband. Therefore, the FrFT of that same signal, with $\alpha \neq 0$, sweeps between $f_{start} = -f_{limit}$ and $f_{stop} = f_{limit}$.

I was wondering, how do I generate a chirp between those two frequencies, $-f_{limit}$ and $f_{limit}$? And, is there any documentation for this?

The one element missing from your problem is the time duration of the pulse - I'll assume it is $T$. The general form of a chirp is
$$
p(t)=\exp\left(j2\pi(f_0t+\tfrac{1}{2}\alpha t^2)\right),
$$
where $f_0$ is the frequency at the start of the chirp and $\alpha$ is the chirp rate. For you case $f_0 =-f_{limit}$. To find $\alpha$ we set the end frequency to be $f_{limit}$ and solve
$$ f_0+\alpha T =f_{limit}$$ which gives:
$$ \alpha = \frac{2f_{limit}}{T}$$

Thus the pulse is given by:
$$
p(t) =\exp \left(j2\pi\left(-f_{limit}t+\frac{f_{limit}}{T}t^2\right)\right)
$$


Effect of aliasing when sampling at higher rate Considering a signal $f(t)=\sin(2\pi t)$ with frequency $10\textrm{ Hz}$:

Middle and bottom plots show the results of sampling at different rates.

I'm trying to understand the case shown at the bottom.

Is the sampling rate too high that it's capturing noise?
What's the name of this effect? May I call it a kind of aliasing caused by high sampling rate?

The name of this effect is Spectral Leakage. Remember the relation
$$\frac{k}{N} = \frac{F}{F_S}$$
where $k$ is the bin number, $N$ is the FFT size, $F$ is the continuous frequency in Hz and $F_S$ is the sample rate in Hz. It can be seen that $k$ varies from $-N/2$ to $N/2-1$ (or from $0$ to $N-1$). So there are only $N$ continuous frequencies $F$ for which you will get a single impulse in discrete frequency domain.

So in the first figure, the frequency of the signal contains an exact integer number of periods within the FFT duration and hence the correlation of other DFT sinusoids with the signal is zero. In the second figure, the signal span is not equal to the integer number of periods that violates the orthogonality requirement. Hence, the correlation of different DFT sinusoids appears as leakage around the true frequency.


Hybrid Sampling Scheme Idea? I am facing a problem which I could not find any robust solution for . Assume we have a signal, $x$, composed of sum of a few sinusoidal.  e.g. $$x=A_1\sin(\omega_1t+\phi_1)+A_2\sin(\omega_2t+\phi_2)+A_3\sin(\omega_3t+\phi_3)+\ldots +A_n\sin(\omega_nt+\phi_n)$$
We have sampled the signal,  ($30\%$) of samples are high resolution 10-Bits and the other $70\%$ are low resolution 3-4 Bits samples.

Also,assume a sub-set of frequencies are known(for example $\omega_1, \omega_5$ and  $\omega_{10}$ out of $\omega_1,\ldots , \omega_5$ are known).

How can we reconstruct the original signal?
Is there any method which can reconstruct the signal from these samples?

Moreover, we have these assumptions:

The sinusoidals have unknown amplitudes ($A_n$)
of course, we know when the high resolution samples occur
We know maximum number of possible independent frequencies in the
summation
and spacing of sampling of high resolution sampling is reconfiguration
and can be programmed to be regular or random (whichever is
suitable?)

For examples, the below is a signal sampled uniformly at Nyquist rate($50\textrm{ Hz}$), with two components of $25\textrm{ Hz}$ and $12\textrm{ Hz}$. Green samples are high resolution samples, red samples are samples with considerable error (low resolution) and the blue curve is the original signal.

Following code provides an example of mentioned data:

Fs=1000; %sampling frequency
t=(0:800).*(1/Fs);
x=0.5*sin(2*pi*300*t)+0.6*cos(2*pi*460*t)+0.2*cos(2*pi*60*t)+0.8*cos(2*pi*40*t)+1.6*cos(2*pi*12*t); %Original Signal
NoisySample=x+0.1*randn(size(x));           %Uniform Noisy Samples
ClearSamples=x(1:10:end); %Clear Samples

I would appreciate any hint and direction.

i will reframe the question a little:

we have

$$\begin{align}
 x(t) &= \sum\limits_{n=1}^N r_n \cos(\omega_n t + \phi_n) \ + \ \epsilon(t) \\
\\
  &= \sum\limits_{n=1}^N a_n \cos(\omega_n t) + b_n \sin(\omega_n t) \ + \ \epsilon(t) \\
\end{align}$$

where

$$\begin{align}
 a_n &= r_n \cos(\phi_n) \\
 b_n &= -r_n \sin(\phi_n) \\
\\
 r_n &= \sqrt{a_n^2 + b_n^2} \\
\\
 \phi_n &= \arg\{a_n - j b_n \} \\
        &= -\operatorname{atan2}(b_n,a_n) \\
\end{align}$$

and $\epsilon(t)$ is the unknown error signal.

first you have to assume that $\epsilon(t)$ contains no components having frequencies of any of the $\omega_n$ in it, otherwise that component will just team up with the component having amplitude $r_n$.

if you know all of the $\omega_n$ in advance, you can calculate all of the $a_n, b_n$ by using an inner product:

$$\begin{align}
a_n &= \frac{\omega_n}{M\pi} \int\limits_{t_0}^{t_0+\tfrac{2M\pi}{\omega_n}} x(t) \cos(\omega_n t) \ dt \\
\\
b_n &= \frac{\omega_n}{M\pi} \int\limits_{t_0}^{t_0+\tfrac{2M\pi}{\omega_n}} x(t) \sin(\omega_n t) \ dt \\
\end{align}$$

$t_0$ can be the most convenient time for you and $M$ should be a large positive integer (and need not be the same $M$ for each frequency component $\omega_n$), but the wider your integration is, the better error rejection you will have.

from the $a_n$ and $b_n$ coefficients, you can calculate $r_n$ and $\phi_n$.


The right way to approach z transform? I am a student learning dsp. I like the subject. I could understand the discrete time signals. When I move into z transform. I could not understand it.

Z transform is the mapping from discrete signal to zplane to make frequency analysis easy. It would be great if users give some tips that could help me to clarify various views on z-tranform on discrete signals.

Text is here to explain about broad answers on z-tranform. I am just looking forward in some tips that can help me understand and approach z-transform in text. Just some lines about what you think about z-tranform will surely help. Thank you mate.

Discrete-time signals are generally written in the time domain as $x[n]$ where $n$ is an integer.

We want to be able to see how discrete-time systems can be used to modify these signals:
$$
y[n] = \sum_{m=1}^M a_k y[n-m] + \sum_{p=0}^P b_p x[n-p]
$$

The above equation is a linear, constant coefficient difference equation.

Such equations can be hard to manipulate for a specific $x$ to figure out what $y$ is.

The $z$ operator is used to attempt to simplify the analysis of such systems. $z$ is the forward-shift operator: it shifts a signal forward by one time instant. $z^{-1}$ is the backward-shift operator: it shifts a signal backward in time by one time instant.

Using this operator, we can rewrite the equation as:
$$
Y(z) = Y(z) \sum_{m=1}^M a_k z^{-m}  + X(z) \sum_{p=0}^P b_p z^{-p}
$$
which, upon rearrangement, we can get:

$$
Y(z)/X(z) = \frac{\displaystyle\sum_{p=0}^P b_p z^{-p}}{\displaystyle1 -  \sum_{m=1}^M a_k z^{-m}}
$$

which is now a discrete-time system that we can analyze in the $z$ domain, rather than the time domain.

One thing this buys us is the ability to figure out what happens when many discrete-time systems are put together: we just multiply the discrete-time system descriptions (transfer functions).  Doing the same thing in the time domain is doable, but very tedious.


What is the real difference between DSP and AI/data science? I am very much into DSP, but lately everybody is talking about "Data science" and "deep learning", and from what I understand the whole thing about data science is to take a huge row of data, and extract meaningful data only from it.

Isn't it what you do with DSP, when you have a signal with noise and you extract the data from a raw of sampling? Or for example when you auto-correlate to find matches between signals?
Do you use the same statistical tools in both cases?

The purpose of both DSP and machine learning is to transform the input signal / data set into more meaningful information. This could be,

removal of noise from an audio recording
location of faces in an image
classification of objects in an image
etc

The main difference as I see it, is that in DSP the transform is designed by the engineer. The engineer will choose a set of signal processing operations that give the desired output. That choice is guided by experience and validated from the results of experiments.

In machine learning the transform is learnt. Typically, this requires a set of training signals with known outputs, over which the system will optimise its transform. The engineer still has to choose a machine learning architecture that can replicate the transform that is needed.

Often the two will be combined. And DSP is kind of machine learning anyway - the machine learning bit is done by the human and the test set is the sum of their existence.


Power Spectral Density and Fast Fourier Transform basic question I am new to signal processing and I am trying to get the basic concepts right.
I have a signal and I performed the FFT of the signal in MATLAB. I also did the PSD of the signal in MATLAB using the psd command.

So I am confused about:

1) When we get the FFT of the signal, It's basically showing what frequency components are inside the signal. So it's a plot of frequency against what? It is amplitude in volt right(took the abs value of the signal)? What if I want to plot dbV against frequency should I just do a mag2dB conversion on the abs value? Is the following right?

Code I tried:

Fs = 16000;
t = 0: 1/Fs : 2*(Fs-1)/Fs;
y1 = sin(2*pi*500*t );
fshift = (fftshift(fft(y1)));
f = -8000 : 1/2 : 7999;
figure;
plot(f,mag2db(abs(fshift)));

2)Now for the PSD, I used the dspdata.psd of MATLAB and I get a plot. On the plot, I get my frequency component but the max amplitude (in dB) is not the same as the max amplitude (in dB) of the FFT plot. Why is that? Is it because dsp is using 10log 10 and mag2dB is using 20 log 10?

3) In the PSD plot I am getting two plots one blue one which is quite visible and a green one, smaller in amplitude, what is this green one?

Code below:

Fs = 16000;
t = 0: 1/Fs : 2*(Fs-1)/Fs;
y1 = sin(2*pi*500*t );
psd1=spectrum(y1,1024);
hpsd = dspdata.psd(psd1,'Fs',Fs); % Create a PSD data object.
figure;
plot(hpsd); % Plot the PSD

The second trace is probably the phase. You can express the spectral density either as a sine term and a cosine term (in phase and out of phase terms) or as a magnitude and phase angle.  In your first example, for your plot you have converted the sin/cos terms into a magnitude, and discarded the phase information.

Yes, log power differs from log voltage by a factor of two.  $ 2 \log{(V)} = \log{(V^2)} $. But different ways of calculating FT and inverse FT, or PSD and inverse PSD, put a constant value of $ \pi $ in different places: If you use different methods of calculation, you sometimes find that you have choosen to normalize the FT differently.


How does harmonic distortion occur in ADCs? We are currently testing an ADC at my workplace and the Total Harmonic Distortion is an important spec. It got me to thinking as to how it occurs in the first place

1) Can there be harmonic distortion if the input signal is a pure sine wave? In other words if the input is theatrically a clean sine wave , can the adc artifically introduce harmonics?

2) ive read quantisation noise can introduce or amplify harmonics of a signal. How does this happen?

Thanks!

There are many distortion mechanisms in ADC but the most important ones are:

"static" distortion due to INL
Track-and-hold distortion

The first one is caused by the fact that the voltage-in-code-out relationship is not perfectly linear (INL/DNL).

An approximation of harmonic distortion based on DNL is given by (Lee, JSSCC 12/2017):

$$
\mathrm{HD3} = -20\log\left( \frac{4}{3\sqrt{3}} \frac{2^B}{\mathrm{INL}_{\mathrm{max}}} \right)
$$

The second one is caused by track nonlinearity. The sampling switch is for example a weakly nonlinear CMOS switch. Its on-resistance is given by (using the simple square law model):

$$
R_{on} = \frac{1}{\mu C_{ox} \frac{W}{L} (V_{DD}-V_{in}-V_t)}
$$

From which is can be seen that the resistance depends on the input signal (see for example Razavi, Data Conversion System Design, page 16).

This distortion can be approximated as

$$
\mathrm{HD3} = 20\log\left(\frac{1}{4} \left(\frac{A}{V_{GS}-V_t}\right)^2 2\pi f_{in} R_{on} C \right)
$$

For a simple CMOS switch this is by far the dominant distortion source and is practically overcome by using a bootstrap sampling switch. However, also this technique has limitations.


signal reconstruction I'm new to DSP so the question might seem a little bit stupid: When using a low pass filter (say Butterworth filter) with Python, I came across an algorithm with which I can reconstruct my signal after filtering (from scipy.lfilter). I'm just wondering how the signal is reconstructed. I thought they used Fourier transform but it's not mentioned.
Thank you!

I agree with A_A. Just take deeper in Butterworth filter design and this article may help you to understand better.


Frequency shifting property Usually in radio systems in order to move the frequency response of the message signal to a new frequency band centered at f, the message signal is mixed with $\cos(2\pi ft)$.

What would change in frequency if we multiply the message signal to  $\sin(2\pi ft)$ instead of  $\cos(2\pi ft)$?

It does the same thing (i.e. shift the frequency to another band), though mathematically it differs a bit. Fourier transforms of the two choices are as shown below:

$$
\begin{array}{|c|c|}\hline
f(x)\cos(ax) & \displaystyle\frac{\hat{f}\left(\xi-\frac{a}{2\pi}\right)+\hat{f}\left(\xi+\frac{a}{2\pi}\right)}{2}\\
\hline
f(x)\sin(ax) & \displaystyle\frac{\hat{f}\left(\xi-\frac{a}{2\pi}\right)-\hat{f}\left(\xi+\frac{a}{2\pi}\right)}{2i}\\\hline
\end{array}
$$

The negative frequencies not only get shifted, but also flipped in polarity

Now for the derivation of these formulas:
Using the Great Mr. Euler's identity: $e^{iax} = \cos(ax) + i \sin(ax)$

$$\cos(ax) = \frac{e^{iax} + e^{-iax}}{2}$$ and
$$\sin(ax) = \frac{e^{iax} - e^{-iax}}{2i}$$

\begin{align}
\textrm{F.T.}\left\{f(x) \cos(ax)\right\} &= \int\limits_{-\infty}^{\infty} f(x) \cos(ax) e^{-i2\pi\xi x} dx\\
&= \int\limits_{-\infty}^{\infty} f(x) \frac{e^{iax} + e^{-iax}}{2} e^{-i2\pi\xi x} dx\\
&= \frac{\int\limits_{-\infty}^{\infty} f(x) e^{-i2\pi(\xi + \frac{a}{2\pi} )x} dx + \int\limits_{-\infty}^{\infty} f(x) e^{-i2\pi(\xi - \frac{a}{2\pi})x} dx}{2}\\
&= \frac{\hat{f}(\xi + \frac{a}{2\pi} ) + \hat{f}(\xi - \frac{a}{2\pi} )}{2}
\end{align}

You can derive the formula for $\sin$.


Sequential processing of uncorrelated measurements in Kalman Filter I'm starting to brush up on the Kalman Filtering I learned a couple decades ago. From what I remember, if you have a measurement vector

$$
z=H x + v
$$
 and the $n$ components of the measurement noise $v$ are uncorrelated with each other, then you can process the measurements sequentially in the KF equations. Is there a good reference somewhere (preferably online) that discusses why this is so? Intuitively it is understandable if each measurement $z$ depends on only one state element in $x$, but it isn't obvious if each measurement of $z$ depends on several elements of $x$.

I think it is evident in both cases. It follows directly from
$$ P(A|B,C,D) = P(((A|B)|C)|D)$$
.i.e.
$$P(x[k]\!=\!\hat{x}[k]|\!\hat{z}_1[k]\!=\!z_1[k],\!\hat{z}_2\![k]\!=\!z_2[k])\!=\! P\Bigg(
  \Big(
      x[k]\!=\!\hat{x}[k]|\!\hat{z}_1[k]\!=\!z_1[k]
  \Big)|\!\hat{z}_2[k]\!=\!z_2[k]
\Bigg)$$
As long as you don't do a prediction update it should be fine.
For nonlinear system, you won't get the same exact estimate, because you will end up with a new dynamics, input, measurement and feedback matrices on every "single-measurement" update


How can noise be added to a signal to obtain a specific SNR? How can I add noise to a recording to obtain a SNR value provided by the user? I am tryin to test the effects of noise on my algorithm. Is there any function or app for doing that kind of task?

Thank you.

Easiest way is to scale energy of noise. Nevertheless you must understand that for some files it is not that simple. Let's say that you have the recording of crying baby in silence. What is then your point of reference? Are you going to add white noise based on ratio between it's energy and average energy across the whole recording (silence+baby)? Or maybe you should calculate the energy in sliding window and then use only sound of interest as point of reference? That is why you should not add noise based on energy of your whole file.

When it comes to some particular examples, then I suggest you to take a look here for mathematical derivation. If you are looking for functions, then please check out the MATLAB awgn function. Nevertheless if you want something even more sophisticated, then you might want to check the Audio Degradation Toolbox. You can find a dedicated routine for adding various types of noise with a given SNR.


Frequency Spectrum of an image I have been studying about Frequency Domain from many days but still i am not able to clear one of the confusion. We say that in frequency domain an image is represented as waves. What I am able to think is that if we take out an intensity profile column wise (i.e one row of intensity values), we have pixel location and at each pixel location we have some intensity value. So using these values with instensity values on y axis and location on x axis, we can draw a wave for one particular row. Now we find DFT coefficients. These DFT coefficients are multiplied with various waves of different frequencies and added. Ultimately we get the original image. Similarly we do it for 2nd row , 3rd row and so on. But when I saw the following image I got confused.

If we take its fourier transform and see its image in frequency domain we get 3 dots in the middle at frequency location 4. I know we get these dots at frequency location 4 because the image has frequency 4. Now the question is why do I get only one row of 3 dots? Why not each row contains 3 dots as I have discussed earlier? Please explain. Thank You.

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


Papers for Teaching an Introductory Course in Image & Signal Processing I'm preparing a digital signal/image processing introductory course and, besides the theory and practices I want the students to read and discuss a couple of research papers during the last one or two classes. It doesn't matter if the papers are new I just want to expose the students some first hand research material in its original format. Which papers could you recommend me? Thanks

I would let them read the paper about the Non Local Means Filter:

Antoni Buades, Bartomeu Coll, Jean Michel Morel - On image Denoising Methods.

The paper is readable and it is a great introductory to the Denoising operation in the context of Image Processing.
Also the Non Local Means is a very decent method (Result wise) even in our time.


Notes or books to quickly get started on the basics of signal processing or wireless communication I'll soon start working in signal processing, specially wireless communication systems. Hence I'd like to get a quick understanding of MIMO and related systems, for example, the topics like:

1) where and how the model $y=Hx + n$ arise,

2) fading, Rayleigh fading,

3) channel capacity, channel state information,

4) other basics one must know

etc.

Note that, I'm not asking for the mathematically rigirous literature; I'm aware I'll need to work with large dimensional random matrices and I'm familiar with the random matrix side of things (I'm a mathematcian), but I'd like to know something bit more basic and applicartion-oriented, for example, where those random matrices actually arise, how the models are formulated etc.

An intrductory (easy to understand) material that'd help me gain all these necessary basics in a week (or say ten hours) would be greatly appreciated. Thank you in advance!

One of the most readable (imho) authors on the communication systems is Simon Haykin and the following book from him would probably address most of the issues you would encounter in a wireless communication system analysis.

Modern Wireless Communication Systems

Similar known authors do have related books, but I guess most compact academic begining would be this one.

It's not a practical book by any comparison, and therefore does not discuss any bit of a practical application of any of the mathematical analysis it intruduces...

If this book contains advanced materials then you may consider reading its pre-requisites from the same author:

Introduction to Analog and Digital Communications

Furthermore you can enhance your Matlab applications on the fundamental communication systems by the following book from the author Proakis.

Contemporary Communication Systems using Matlab


How to design an all pass type filter with inverted phase response? I want to model a "negative dispersion" in waveguides.

Is it possible to design an all pass filter with inverted phase response? (e.g. -pi at 0 and pi at Nyquist)?

I feel :) this could be possible in complex domain (but don't know how). I could use a Hilbert transformer to complexity real signal at first.

Sort of. The inverse of an allpass filter is also an allpass filter and can simply be calculated by time flipping the impulse response. However that makes the filter non-causal. So you need to add bulk delay to make it causal again. As long as you don't need to do real time processing, that's typically not a problem.

EDIT: add a specific example. The key is to do the filtering in the time domain. You need to truncate the IIR response to an FIR response using "enough" samples.

% create an arbitrary allpass filter
b = [-.5 1]; a = [1 -.5];

% calculate the impulse response
n = 128; d0 = zeros(n,1); d0(1)= 1;
h1 = filter(b,a,d0);

% time flip and shift by 1 to keep the first sample in the same spot
h2 = circshift(flipud(h1),1);

% plot frequency responses
fh1 = fft(h1);
fh2 = fft(h2);
phases = unwrap([angle(fh1(1:n/2+1,:)) angle(fh2(1:n/2+1,:))]);
plot(phases);


Recommended device for testing my own SDR? I'm trying to implement my own SDR using Matlab. I was wondering what kind of device I could buy in order to test it. I know I can use my own soundcard but I was thinking of EM waves, not sound. Is there a device capable of functioning as transmitter and receiver and be compatible with Matlab? Just wondering, I have found very cheaps ones (like RTL-SDR) but they only work as receivers and I'm not sure if they're Matlab-compatible. Thanks!

For both receive and transmit SDR experimentation/testing at RF/HF/VHF/etc. frequencies, I use a LimeSDR Mini (far less expensive than an Ettus). But I stream from/to it using C programs and files, not direct from Matlab.  (Is it possible to write a C stub to pass data back and forth with Matlab?)

There is also the ADALM-Pluto kit from Analog Devices, perhaps suitably priced for educational environments.  There seems to be some sort of Matlab support for it.

Some amateur radio transceivers (Icom, et.al.) have stereo soundcard (audio IQ) ports for conversion from/to (licensed) amateur radio frequency band RF.

Added: Another solution for amateur radio HF band SDR experimentation is the Hermes Lite 2 open source kit transceiver


How to filter my noisy square wave so that i can recover my signal using the Chebyshev Filter I have a noisy square signal as input but it has a lot of noise.
I should use a Bandpass filter to recover my signal.
I know that the Chebyshev Filter is a bandpass filter; but it doesn't work. the data is still noisy after filtering.

Chebyshev Filter:
The Chebyshev filter  gives a sharper cutoff than a Butterworth filter in the pass band. A disadvantage of the Chebyshev filter is the exterior of gain minima and maxima below the cutoff frequency. Due to the ripples in the pass-band, it is not used in the audio systems. Though it is far better in some applications where there is only one frequency available in the pass band, but numerous other frequencies are required to eliminate.

 % I use A picoscope2000 serie for data acquisition. (MATLAB)
 ps2000_getdata;

 % The data is saved a cha_a (MATLAB)
 cha_a = (bufferChA/1000);

 % Here i use the Chebyshev bandpass filter to filter the noisy data(MATLAB)
 [A,B,C,D] = cheby2(10,40,[900 1100]/1500);
 d = designfilt('bandpassiir','FilterOrder',20, ...
'StopbandFrequency1',900,'StopbandFrequency2',1100, ...
'StopbandAttenuation',40,'SampleRate',8000);

  y= filter(d,cha_a); % convolution

  sos = ss2sos(A,B,C,D);
  fvt = fvtool(d,y,'Fs',8000);
  legend(fvt,'cheby2','designfilt')

This is the result of the magnitude response :

And this is the impulse response :

My input data is suppose to be a squared data.
And when you look at the impulse response of the output, it does not look like the impulse response of a square signal. There is stil noise. How can i adjust my filter ? Or did i used the filter well ?

Bandpass filtering won't do what you're after.  If I take a noiseless square wave with the spectrum in red, and design a bandpass filter for it (response in blue):

and then apply this filter to the original (light blue) square wave, then I get the dark blue response.

I suspect that a better approach might be to use a median filter.

If I apply some gaussian noise to the square wave (blue line), and then apply a median filter (orange plot) then the result is below.

These plots were generated in python; I don't have easy access to Matlab. Code below.

from scipy import signal
import math
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0,20*math.pi,1000)
sq = scipy.signal.square(t)

b, a = signal.cheby1(4, 0.1, [0.01, 0.2], 'bandpass')

plt.rcParams['figure.figsize'] = [20, 10]
plt.figure(1)
w, h = signal.freqz(b, a)
w2, h2 = signal.freqz(sq, 1)
plt.plot(w, abs(h)*600)
plt.plot(w2, abs(h2), color='red')

plt.figure(2)
plt.rcParams['figure.figsize'] = [20, 10]
plt.plot(t,sq)
sq_f = scipy.signal.lfilter(b,a,sq)
plt.plot(t,sq_f,color='blue')

sq_noisy = sq + np.random.normal(0,1,1000)
plt.figure(3)
plt.plot(t,sq_noisy)
sq_noisy_f = scipy.signal.medfilt(sq_noisy,21)
plt.plot(t,sq_noisy_f,color='orange')


IMU Speed Tracking Through Known Path I am new to signal processing and Kalman Filtering here. Thanks for your help.

I working with an IMU for a tracking project where the IMU moves throw a known path but at an unknown speed (within limits), the objective being tracking the speed at which the IMU moved through the path. To my surprise, I couldn't find anything online on IMU tracking through a predetermined path.

The idea of using the known, 'predetermined trajectory', of the IMU as a constrain makes sense to me, using the path that we know that the IMU is going to travel through as a form of container. But I don't know if there is a flaw in my reasoning or if I'm missing something.

I thought a Kalman filter could be used to merge the accelerometer and gyroscope data and have it constrained to the trajectory as a way to limit sensor noise, but I rly can't get my head around how to do this! Maybe this is an optimization problem where the measured 3D trajectory is optimized into the 'predetermined trajectory' instead?

Ideally this could be applied to any kind of movement to which the trajectory is known, so I'm imagining this 'predetermined trajectory' as a 3D curve.

You should parameterize the path as a parameter of time.
You can do that off line with accurate measurements of the path.
Then use Non Linear Least Squares to find the best match between the reads of the IMU to the model.

From there you'd be able to match the speed.


How to calculate a delay between two signals in frequency domain? Calculation of a delay between two signals is generally done by finding the maximum in the cross-correlation of the two signals (in time domain):
$$
\tau_{delay} = \arg\max_\limits{t}(f(t) \star g(t))
$$
A cross-correlation can be easily computed in the frequency domain by pointwise multiplication of the Fourier transform of one signal with the complex conjugate of the Fourier transform of the other signal:
$$
f(t) \star g(t) = f^*(-t)*g(t)=\mathcal{F}\{F\}^*\cdot \mathcal{F}\{G\}
$$

Now I wonder how I can compute the cross-correlation peak in the frequency domain, thus without the need of an inverse transform. I can imagine two easy cases:

The case for two sinusoidal signals seems trivial, as the (smallest) delay then agrees with the phase (divided by the angular frequency) of the Fourier component (of this pointwise multiplication) for the respective frequency.
This can be generalized to two identic signals with a fixed delay, as the (unwrapped) phase response equals then the angular frequency multiplied by the delay, for all Fourier components. To avoid ambiguities, let's assume that the delay is smaller than half the period of the highest frequency that is considered.

However, when the two signals are not exactly identic, I suppose some kind of "weighed sum/integration" should be taken of the phase response, where the weighting factor scales with the relative contribution of the respective frequency, i.e., with the magnitude of the respective Fourier component. However, I am not able to find the exact formula rigorously.

Some "internal brainstorming" that might lead to the solution: if we consider one of the two functions ($f$) as a "basis" function, and project the other signal ($g$) on this function, the result of this projection (inner product) should contain the delay we are looking for. I guess this is only possible if the "basis" is first constructed as the analytic representation of $f$, such that phase differences can be taken into account. I also think that this inner product should be conserved after Fourier transform.

Calculating a delay between two signals is not as simple as you described and it cannot be generalized as easily either.

  The case for two sinusoidal signals seems trivial, as the delay then agrees with the phase (divided by the angular frequency) of the Fourier component (of this pointwise multiplication) for the respective frequency.

This is not true in the general case. Calculating the delay from a single sinusoidal component is inherently ambiguous because of the signal's periodic nature. For a single phase $\phi$, there are infinitely many valid delays $\tau$.

$$
\tau = \frac{2\pi \cdot n + \phi}{2\pi f}, n \in \mathbb{Z}
$$

  This can be generalized to two identic signals with a fixed delay, as the phase response equals then the angular frequency multiplied by the delay, for all Fourier components.

This isn't true either. As noted above, the phase response plus a multiple of 2π equals the angular frequency multiplied by the delay, for all Fourier components. Even worse, this is a different multiple for each Fourier component in the general case.

You want to find the delay $\tau$ that fulfills this equation for all Fourier components, or at least for most of them, if the signals are not exacly identical.

The inverse fourier transform solves this problem in a brute-force manner. For each candidate delay, it calculates how strong the measured phase response agrees with it. Just pick the maximum and you're done. There is no better algorithm for the general case because the problem is inherently ambiguous and many (or even all) solutions might agree equally strong. Obviously, there cannot be a method to calculate the solution of such a problem directly.

Things change if you want to calculate the delay with higher precision than one sample. You could just zero-pad the phase response before calculating the inverse fourier transform, but this is computationally inefficient. A better approach is to calculate a coarse solution first, using an unpadded inverse fourier transform and then refine it using Newton's algorithm. Newton's algorithm is perfecly fine here because

it never requires the calculation of a delay from a phase response (which would be ambiguous). You can get away with calculating agreement scores for a given phase response and a candidate delay.
convergence to a local maximum is perfectly fine if all you wanted to do is to refine the global maximum.

There is no non-brute-force algorithm to find the global maximum of an arbitrary function. You are trying to find the global maximum of the cross correlation function, so unless you can guarantee certain properties of the signals (such as convexity or the like) or you're happy with a local maximum, there will be no way around an inverse fourier transform.


What is the true meaning of a minimum phase system? What is the true meaning of a minimum phase system? Reading the Wikipedia article and Oppenheim is some help, in that, we understand that for an LTI system, minimum phase means the inverse is causal and stable. (So that means zeros and poles are inside the unit circle), but what does "phase" and "minimum" have to do with it? Can we tell a system is minimum phase by looking at the phase response of the DFT somehow?

The relation of "minimum" to "phase" in a minimum phase system or filter can be seen if you plot the unwrapped phase against frequency.  You can use a pole zero diagram of the system response to help do a incremental graphical plot of the frequency response and phase angle.  This method helps in doing a phase plot without phase wrapping discontinuities.

Put all the zeros inside the unit circle (or in left half plane in the continuous-time case), where all the poles have to be as well for system stability.  Add up the angles from all the poles, and the negative of the angles from all the zeros, to calculate total phase to a point on the unit circle, as that frequency response reference point moves around the unit circle.  Plot phase vs. frequency.  Now compare this plot with a similar plot for a pole-zero diagram with any of the zeros swapped outside the unit circle (non-minimum phase).  The overall average slope of the line with all the zeros inside will be lower than the average slope of any other line representing the same LTI system response (e.g. with a zero reflected outside the unit circle).  This is because the "wind ups" in phase angle are all mostly cancelled by the "wind downs" in phase angle only when both the poles and zeros are on the same side of the unit circle line.  Otherwise, for each zero outside, there will be an extra "wind up" of increasing phase angle that will remain mostly uncancelled as the plot reference point "winds" around the unit circle from 0 to PI.  (...or up the vertical axis in the continuous-time case.)

This arrangement, all the zeros inside the unit circle, thus corresponds to the minimum total increase in phase, which corresponds to minimum average total phase delay, which corresponds to maximum compactness in time, for any given (stable) set of poles and zeros with the exact same frequency magnitude response.  Thus the relationship between "minimum" and "phase" for this particular arrangement of poles and zeros.

Also see my old word picture with strange crank handles in the ancient usenet comp.dsp archives: https://groups.google.com/d/msg/comp.dsp/ulAX0_Tn65c/Fgqph7gqd3kJ


Getting introduced to signal processing for remote sensing I work mainly in GIS (geographical information systems) but I would like to deepen my knowledge of remote sensing (analyzing imagery acquired remotely from airborne or satellite sensors). I was advised to learn as much about signal analysis as possible. I am assuming that signal processing and analysis refer to the same techniques.

Because signal processing seems to be applied to all sorts of different areas, I'm not quite sure where to start. I should specify that beyond two university calculus courses (on which did quite well), I do not have a background in physics or math.

Could anyone help me plan a sort of "road map" to learning signal processing, preferably applied to remote sensing? Thank you.

@AndrewP. You have a long road ahead you. But the road is well lit and the rewards will are greater and greater as you travel along. You need to learn digital signal processing (DSP). May I humbly suggest you buy the two books "The Scientist & Engineer's Guide to Digital Signal Processing"

http://www.amazon.com/Scientist-Engineers-Digital-Signal-Processing/dp/0966017633/ref=sr_1_2?ie=UTF8&qid=1438652461&sr=8-2&keywords=scientist%27s+guide+signal+smith

and "Understanding Digital Signal Processing."

http://www.amazon.com/Understanding-Digital-Signal-Processing-3rd/dp/0137027419/ref=ntt_at_ep_dpi_1

They're widely considered to be the best ways to learn DSP for beginners. Do NOT sign up for any free Internet online DSP courses!  They will only discourage you at this point in your learning. You'll need to obtain either (1) a 'Student' version of Matlab software or (2) Octave software so you can experiment with DSP on your computer. Do not take any math classes at this point--you can learn selected mathematical topics as necessary as you travel your long road. Andrew if you have enough energy and enthusiasm, to quote Susan B. Anthony, "Failure is impossible."


Determining the autocorrelation sequence from an AR model I have the following equation:

$$x(n)=\frac{14}{24}x(n-1)+\frac{9}{24}x(n-2)-\frac{1}{24}x(n-3)+w(n)$$
where, $w(n)$ is a stationary white noise process with variance $\sigma^2_w$

Now, I want to determine the auto correlation sequence $\gamma_{xx}(m)$, for $m = 0,1,\ldots,5$

Now, I know that the autocorrelation function for AR process which is given below:

$$\gamma_{xx}(m)= \begin{cases}\displaystyle-\sum_{k=1}^p a_k\gamma_{xx}(m-k), & m > 0\\
\displaystyle  -\sum_{k=1}^p a_k\gamma_{xx}(m-k) + \sigma_w^2, & m = 0\\
\gamma_{xx}^*(-m), & m < 0\\
\end{cases}$$

Here, I already found the coefficients of ${a_k}$ which are $a_1=-\frac{14}{24}, a_2=-\frac{9}{24}, a_3=\frac{1}{24}$.

Now, I don't get, how can I calculate $\gamma_{xx}(0)$ to $\gamma_{xx}(5)$ recursively?

Please someone tell me how can I solve  $\gamma_{xx}(0)$ and $\gamma_{xx}(1)$? Then I can do the the rest.

Looking at your definition of $\gamma_{xx}(m)$, it seems that you have 4 equations with 4 unknown (for $m=0,\ldots,3$):

\begin{align}
24\gamma_{xx}(0) &= 14\gamma_{xx}(1) + 9\gamma_{xx}(2) - \gamma_{xx}(3) + 24\sigma_w^2  \\
24\gamma_{xx}(1) &= 14\gamma_{xx}(0) + 9\gamma_{xx}(1) -\gamma_{xx}(2) \\
24\gamma_{xx}(2) &= 14\gamma_{xx}(1) + 9\gamma_{xx}(0) - \gamma_{xx}(1)\\
24\gamma_{xx}(3) &= 14\gamma_{xx}(2) + 9\gamma_{xx}(1) - \gamma_{xx}(0)
\end{align}

Solving the above system of equations is straight forward, and after that you can easily obtain $\gamma_{xx}(m)$ for $m=4,5$


Fourier transform in Matlab and hermitian symmetry According to the conjugate symmetry property of Fourier transform, shouldn't the following command not return 1 (=true):

x=imread('cameraman.tif');
ishermitian(fft2(x))

However it does not (returns 0).

Is it due to some rounding error? Something else? Otherwise, how would you check that it is (avoiding to write a painful code with for loops etc)?

In the 2D (continuous) Fourier domain, the Hermitian symmetry (for functions) writes (conjugate symmetric with respect to the origin):

$$F(-u,-v) = \overline{F}(u,v)\,.$$

For matrix $A_{i,j}$, we should have:

$$ a_{i,j} = \overline{a}_{j,i}$$

which means that the diagonal is real. As you can see, Hermitian matrices and Hermitian functions are slightly different. There is something in the diagonal that needs a center. But diagonals of matrices are real.

Plus, the FFT output matrix is not straightforwardly organized like the (continuous) Fourier transform. Friday night here, I am coming back later on FFT symmetries.


What are Simulated B-mode and phantom images? Could someone please explain what simulated B-mode and phantom images are? I've seen them in several IEEE papers. How can one create these images?

What is the formula for signal-to-noise ratio and contrast-to-noise ratio?
Please share MATLAB code for the above ratios.

B-mode is the most elementary way to form an image from the Ultrasound modality. You basically shift the beam of an A-Mode, either mechanically or electronically. To simulate B-mode, you need a three dimensional model of tissue density and a sound propagation model. There are basically two methods to achieve it, you either do ray-tracing (if the wavelength is much smaller than the distance) or you solve the wave equation directly. In either of these cases you simulate ultrasound propagation through the tissue and re-construct the image that a real ultrasound machine would give you.

A "phantom", more generally, is an item that has been fully characterised and its measurements are known. It is used for calibration. So, for example, you know that Water should return a zero in the Haunsfield scale and Air should return -1000. If you scan a phantom and the image does not agree with its specifications, the machine needs to be re-calibrated / repaired.

The simplest example of a phantom in "everyday life" is that gray-scale ribbon that you see in some screen calibration utilities that direct you to tune your monitor until you can see all the bands in the ribbon, or until some patterns appear / disappear.

A software phantom is exactly the same thing but it is used to "calibrate" or assess an algorithm. You are probably talking about the Shepp-Logan phantom which contains a number of fully characterised ellipses.

If you are working in tomographic image reconstruction then you can compare the image that your technique reconstructs against the "real" image and make inferences.

To create these images, all that you have to do is create an image with shapes of known characteristics (as you can see from the Wikipedia page on the Shepp-Logan phantom).

As far as the "...formula for signal-to-noise ration and contrast-to-noise ratio..." is concenred, it can be found on Wikipedia. It really depends on what you are trying to measure and how and yes, you will need a phantom for that. For more information please see this link.

Image quality assessment is a field on its own in medical imaging.

Hope this helps.


Using the Inverse Filter to Correct a Spatially Convolved Image (Deconvolution) As part of a homework assignment, we are implementing the Inverse Filter. Degrade an image then recover with an Inverse Filter.

I convolve the image in the spatial domain with a 5x5 box filter. I FFT the filter, FFT the degraded image, then divide the degraded image by the filter. Inverse FFT the result into an image and I get garbage.

If I FFT the image, FFT the filter, multiply the two, divide that result by the FFT'd filter, obviously I get very close to the original image. ((X*Y)/Y ~== X)

I have an inkling the math is not as simple as "spatially convolved == FFT multiplication".

What is is the correct way to use the Inverse Filter? I have the exact kernel used degrade the image. I'm not adding any noise.

Bovik's textbook, The Essential Guide to Image Processing is almost completely dismissive of the Inverse Filter. Gonzalez&Woods is a bit more hopeful but almost immediately skips to the Wiener Filter.

I have a similar question on stackoverflow.com
https://stackoverflow.com/questions/7930803/inverse-filter-of-spatially-convolved-versus-frequency-convolved-image

(This questions should also be tagged [homework] but the tag doesn't exist yet and I haven't the rep to create it.)

EDIT. For some of the great suggestions below.
@dipan-mehta Before I FFT, I am padding the convolution kernel to the same size as the image. I'm putting the kernel into the upper left. I ifft(ifftshift()) then save to an image and I get a good result. I've done the ifft(ifftshift()) on both the kernel and the image. Good(ish) results. (Images are in my https://stackoverflow.com/questions/7930803/inverse-filter-of-spatially-convolved-versus-frequency-convolved-image question.)

@jason-r is probably correct. I don't understand the mathematics of the underlying convolution + transform. "Deconvolution" was a new word for me. Still have much to learn. Thanks for the help!

My solution for the homework assignment is to do everything in the frequency domain. I spoke with the professor. I was making the assignment harder than necessary. She wanted us to add noise then try the Inverse Filter, Wiener Filter, and Constrained Least Squares Filter. The point of the exercise was to see how the filters handled noise.

There's a couple subquestions that I'll address separately:

Convolution in the spatial domain (or correspondingly in the time domain for time-sampled signals) is equivalent to multiplication in the frequency domain. In sampled systems, there are some subtleties to boundary cases (i.e. when using the DFT, multiplication in the frequency domain actually gives you circular convolution, not linear convolution), but in general, it really is that simple.

Pure inverse filtering is almost never the right solution in practice. In most cases, you don't have access to the exact filter that has been applied to your data, so you can't simply invert it anyway. Even if you do know the filter, then it's still problematic. Consider the fact that the filter may have zeros at certain spatial frequencies; if it does, then after applying the filter to your image, all information at those frequencies is lost. If you naively invert that filter, it will have infinite (or at least very high) gain at those nulls. If you then apply the naive inverse to an image that has any additive content at those frequencies (e.g. noise, which is likely to be the case), then that unwatned component will be greatly amplified. This is generally not desirable.
This inverse-filtering problem is very similar to equalization in communications systems, where this phenomenon is referred to as noise enhancement. In that context, the inverse-filter approach is referred to as a zero-forcing equalizer, which is rarely actually used.

The area that you're exploring is known more generally as deconvolution. As a general rule, deconvolution is a tricky operation. Even if you know the exact filter that was applied and want to undo it, it's not always that easy. As you noted, the inverse filter approach is usually brushed aside in favor of a Wiener filter or some other structure that aims not to exactly invert the system, but instead to estimate what the input to the system was while minimizing some error criterion (minimizing mean-squared error is a common goal). As you might expect, applying a Wiener filter to this problem is referred to as Wiener deconvolution.


iir biquad filter not working I am trying to implement a butterworth bandpass filter
My filter requirement is Fs= 500Hz, Fc1=0.4Hz, Fc2=4Hz, Oreder=2.

Below is the code which is implemented on 16bit-mcu and is called for every new sample
The coefficients are calculated from http://www-users.cs.york.ac.uk/~fisher/mkfilter

static int32_t xv[4+1], yv[4+1];
int32_t my_filter(int32_t sample)
{
    xv[0] = xv[1]; xv[1] = xv[2]; xv[2] = xv[3]; xv[3] = xv[4];
    xv[4] = sample / 1.989695029e+03;
    yv[0] = yv[1]; yv[1] = yv[2]; yv[2] = yv[3]; yv[3] = yv[4];
    yv[4] = (xv[0] + xv[4]) - 2 * xv[2]
                 + ( -0.9380260977 * yv[0]) + (  3.8116138294 * yv[1])
                 + ( -5.8091337665 * yv[2]) + (  3.9355459731 * yv[3]);
    return yv[4];
}

This is the output of code
What i want is output similar to scipy_filter.

Change all data types to "float"
Split the filter into two second order sections
If you can't afford the floating point, you need to learn fixed point math, which will require significant effort


What is the delta parameter in $\tt scal2freq$? I am running complex morlet (cmor1-1.5) wavelet transform. I want to have the transformation for the range of $[1,150]\textrm{ Hz}$ with a step of $1\textrm{ Hz}$.
I want to use scal2freq Matlab function to convert scales to pseudo frequencies but I don`'t know how to define the delta.

The delta argument in the scal2freq MATLAB function is the sampling period that relates to the data sampling frequency (i.e. delta = 1/Fs).  A good description of this function can be found here: http://au.mathworks.com/help/wavelet/ref/scal2frq.html it also includes an example including how to calculate and define the delta argument.  I hope this helps!


Does sum of random variables follows gaussian distribution although each have different mean and variances? I understand Central limit theorem, but cannot understand the result when it applies to diverse situation.
As far as I understand, Central limit theorem states follows:

No matter what random variable has its distribution, if we pick n samples and mean (or just sum) them and do this many time, those values will follow normal distribution (if n is big enough).

And there is a good example for this: a Brownian motion.
If we measure the distance that particles moved suspending in a fluid, distance will follow normal distribution because sum of random hit(movement) of fluid particle will follow normal distribution according to Central limit theorem.
However, for central limit theorem applies, random variables which are merged should follow all the same means and variances. In upper example, all random particle movement will follow same mean and variances, because they have same temperature.
But as far as I understand, there are many cases where noise follows normal distribution but there components are not necessarily follow same mean.
(S_n = (X_1 + X_2 + ... X_n) / n and S_n follows normal distribution but X_i s have all DIFFERENT means)
For example, electrical noise will be composed of many other things which are not necessarily follow equal to each mean. But why electrical noise can be approximated as gaussian noise? Not only electrical noise, why many other things are treated as normal distribution(S_n follows normal distribution) although their components don't have equal mean?(X_i s have all DIFFERENT mean and variances)
Thank you for reading!

You can add random variables with different means and variances. For example, look at two continuous uncorrelated normal-distributed random variables $x_1$ and $x_2$ with means and variances $\mu_1, \sigma_1^2$ and $\mu_2, \sigma_2^3$
Then $x_3 = x_1 + x_2 $ is also a normal-distributed random variable with $\mu_3 = \mu_1 + \mu_2$ and $ \sigma_3^2 = \sigma_1^2 + \sigma_2^2 $. For uncorrelated random variables the means and variances just sum.
The distribution of a sum of two continuous random variables is the convolution of the individual  distributions. That's actually the core of the central limit theorem. If you convolve a Gaussian with a Gaussian you get another Gaussian (see for example https://jeremy9959.net/Math-5800-Spring-2020/notebooks/convolution_of_gaussians.html) . If you keep convolving reasonably looking functions with each other eventually it will approach a Gaussian.
A simple example is  summing a few uniformly distributed random variables. The first one has a rectangular distribution. The second one is triangular. The third one already shows a nice round belly in the middle and once you are up to summing 10 or 12 it looks very much Gaussian. That's in fact one method how gaussian random number generators are implemented.
For the Central Limit Theorem, there is no requirement for random variables to have the same mean or variance.