How to achieve image straightening as cam scanner does? Is it possible to achieve image straightening as cam scanner does

I mean selecting a non-rectangular region and converting it to a rectangular region as an image. I can only do re sizing and crop of rectangular region of image.

Geometric Transformation

Supose that $f$ is an image, defined over $(w,z)$ coordinate system, undergoes geometric distortion to produce an image $g$, defined over $(x,y)$ coordinate system. To perform this operation is used this

$$(x,y) = T\{(w,z)\} $$

So, one of most commmonly used forms of spatial transformation is the affine transform (Wolberg [1990]). It can be written in the matrix form as

$$ [x \ y \ 1] = [w \  z \ 1] \ \mathbf{T} = [w \  z \ 1]
\left[ { \begin{array}{ccc}
         t_{11} & t_{12} & 0 \\
         t_{21} & t_{22} & 0 \\
         t_{31} & t_{32} & 1 \\
        \end{array}} \right]$$

This transformation can scale, rotate, translate or shear a set of points, depending on the values choosen for the elements of $\mathbf{T}$

Your case: Perspective Transformation

For perspective transformation, you need a 3x3 transformation matrix. Straight lines will remain straight even after the transformation. To find this transformation matrix, you need 4 points on the input image and corresponding points on the output image. Among these 4 points, 3 of them should not be collinear.

Implementations

Python

OpenCV

OpenCV has a lot of methods to perform image warping and geometric transformation an you can get the examples in the docs page

scikit-image

scikit-image also have many examples in the documation page.

For python just choose the module more comfortable to you or use both.

MATLAB ® Image Processing Toolbox

As the others examples the Math Work's software have a lot of functions to perform geometric transformation. Beware that is an commercial software and it not have a strong programming language if you want make a little software.

References

Wolberg [1990] G. Wolberg. Digital Image Warping . IEEE Computer Society Press, 1990.

Gonzalez [1987] Gonzalez, R.C., P. Wintz, 1987. Digital Image Processing. 2nd Ed., Addison-Wesley, Reading:MA.

Futher Reading

This lecture is a good introduction to this process


What's the difference between calculating the FFT of a constant sound for 5 seconds versus every 1 second and taking the average? I'm interested in learning the frequency spectrum of a constant sound. However, there is background noise. My plan is to record the sound for 1 second, take the FFT, and then doing this 5 times and taking the average (for each frequency bin). This should decrease the variance of the magnitudes across frequencies.

But I wonder what's the difference (ie, pro/con) of this versus just recording for 5 seconds and calculating the FFT?

Thanks

I would say that you answered part of your own question by noting the reduction of variance with averaging.

Most books tend to present signals that are either random or deterministic but  signals can be somewhat of both character to varying extent.

If your signal tends towards the deterministic, and high SNR, there isn’t a lot variation to reduce so the higher resolution of the longer FFT would probably be more preferable in most applications.

If you have a low SNR tone, you want to do as much averaging as other considerations would permit to reduce that noise.

One can also overresolve with long FFTs.  Pure steady tones are rare. Real tones have some bandwidth and you want your FFT bins to match that bandwidth.

Averaging has some other advantages like reducing the down stream data rates like if you are displaying real time STFT displays.

Depending on your application, there are other more specific tradeoffs between averaging and FFT length.


Beginner Question: Extracting Features for Image Classification I'm a Python programmer, but a beginner to image processing :) Apologies if this isn't a suitable question for the forum, happy to rewrite or move.

I want to create a supervised classifier to distinguish between the two images in this album (or more generally, with more training data, this scene containing some pavement for available parking, and this scene containing no pavement).

My question is about how to extract features, using Pillow or scikit-image, that would be appropriate for classification.

Right now, I have skeleton code for an SVM classifier:

data = []
classes = []
for imagefile in glob.glob('./img/training/*/*.bmp'):
    data.append(extract_features(imagefile))
    if 'nospace' in imagefile:
        classes.append('nospace')
    else:
        classes.append('space')
clf = svm.SVC(gamma=0.001, C=100.)
clf.fit(data, classes)

And I have a very simple feature extraction method, which just returns a colour histogram for the image:

def extract_features(imagefile):
    im = Image.open(imagefile)
    width, height = im.size
    im = im.crop((0, int(height * 0.4), width, int(height * 0.8)))
    # Improve on just returning a colour histogram?
    return im.histogram()

The problem is that this doesn't give great results - particularly if the car is nearly the same colour as the tarmac, as in my example album.

I looked at edge detection. This seems to be quite easy to do with Pillow:

im.filter(ImageFilter.FIND_EDGES)

But just passing this in to the classifier also doesn't give great results - I'm guessing because I should turn it from a matrix of pixels into something more meaningful?

Could anyone suggest a good approach for extracting features for this problem? I think ideally I want to detect "blobs", i.e. cars.

You are just using the histogram as your features, yet your image contains very slight differences. This is not discriminative enough for SVM, you should focus on better features. In this case, you might be better of with subtraction, change detection and etc.


Noise removal in discrete signal How do i remove these kinds of noise from the discrete signal in python?

A simple approach would be to implement a filter using scipy.signal.lfilter, where plenty of documentation exists in python on how to use that function.
To do this effectively you need to first define the signal bandwidth of interest relative to the sampling rate. As long as the sampling rate is sufficiently more than twice the bandwidth, there will be an opportunity for noise reduction. If the noise is truly random, and no knowledge of the signal (equally random) we cannot reduce noise that is occupying the same bandwidth if frequency and duration in time as the signal, but the filter described would be effective in rejecting out-of-band noise.


FFT of resultant of signals If one has two signals (say, two acceleromters mounted perpendicularly) and a piece-wise resultant acceleration signal is determined, it appears that frequency content information cannot be determined on the resultant signal using an FFT.

Can someone comment on why this is so?

It may be that the process of determining the resultant converts the negative content in the individual signals through the process
$$r(i) = \sqrt{x(i)^2+y(i)^2}$$ which "loses" information (not sure if that is the best way to describe it). This same effect can be seen if one takes a single acceleration signal and processes the FFT of the absolute value of that signal.

Is there any way of still being able to assess frequency content of a resultant vector?
Is it even worth pursuing?

First, you should note that the FFT refers to an algorithm for the fast computation of the DFT taking advantage of its symmetry.

The signals you have are:

where $\mathbf{a}_y$ and $\mathbf{a}_x$ can both have negative a positive values depending on the direction along the $x$ and $y$-axis. These signals are the output of your accelerometers.

What you call "a piece-wise resultant acceleration", calculated like
$$r(i) = \sqrt{x(i)^2+y(i)^2},$$ is just $\rVert \mathbf a\lVert$, the magnitude of a vector quantity $\mathbf{a}$, which doesn't show the direction.

From NASA webpage, about vector components:

  Why go to all this trouble? Because, in aerospace, we are often
  dealing with forces and forces are vectors. Breaking a single vector
  force into several components allows us to study the resulting motion
  much more easily.

In conclusion, performing a DFT on your vector quantities $\mathbf{a}_x$ and $\mathbf{a}_y$ will give you frequency information both for about magnitude and the direction together, while performing a DFT on your scalar vector $\rVert \mathbf a\lVert$ will just give you the frequency domain of the magnitude, so yes, you are loosing information here.


map float to fixed point I have  designed a system in simulink and works well through simulation. In order to convert to fixed point, I tried fixed point advisor with little success. The reason being that the coefficients on the float model are larger than 1. I have two of them, a1=0.021 and a2= 3600. My input is from an 16 bit ADC (Q15). So how to scale everything within the range -1 to +1?

Fixed point processing is not for the faint of heart. It typically requires managing scaling factors, clipping point, signal to noise ratios, underflow, limit cycles, rounding behavior at each single point of your processing chain.
Your best shot at dealing with different data ranges is to rigorously use the Q-number format for each data, state and coefficient in your processing (see https://en.wikipedia.org/wiki/Q_(number_format)) and do the correct multiply & shift  & round/truncate at each operation. This may require statistical or worst case analysis to determine the clipping points and best Qs for each variable.
In order of complexity, you can consider

Convert your fixed point data to floating point, run the algorithm in floating point, and then only limit and convert back to fixed point as the last step.
Same as item one but if you can't do float consider using 24-bit or 32-bit fixed point data instead. This creates a extra margin and headroom and makes life A LOT easier
Stick with 16-bit all around. That means you have to fight for every single dB of signal to noise ratio at every step of your processing. Depending on what you are doing you may consider selectively using double precision for critical state variables or operations.


How to find average and power of a signal I need to find the average and the power of this signal:
$$x(n)=\sum_{k=1}^{\infty}2^{-k}e^{j2{\pi}kn}$$
The problem is that the summation starts at 1 and not at 0, and a part of that how can I find the average of a signal expressed in this way?

No solutions to homework problems, but here are two hints:

what is $e^{j2\pi nk}$ for integer $k$ and $n$?
do the values of $x[n]$ depend on $n$?


How is $-6$ dB pressure amplitude is the same as FWHM? I keep reading many literature, and in most cases $-6$ dB is referenced as the full width at half maximum (FWHM). Mathematically, how is that possible from the definition of dB?


What is difference between OFDM and GFDM? I am trying to understand the difference(s) between OFDM (Orthogonal Frequency Division Multiplexing) and GFDM (Generalized Frequency Division Multiplexing) which are used as a multicarrier modulation techniques in wireless communication. As far as I know the OFDM uses the orthogonal subcarriers and the GFDM is with the non-orthogonal subcarriers when designing a transmitter. The GFDM is possibly going to be the modulation technique for 5G systems as it has been discussed in a research area.

Could anyone please explain further about GFDM design and its differences from OFDM?
What is the purpose of using GFDM when designing a transmitter? The subcarriers are created just after the IFFT?

In GFDM, how does the IFFT block work when creating the subcarriers?

What about the subcarrier pulse shaping in GFDM?

Thanks in Advance!

There is a great interactive demonstration of GFDM by Xiaojie Wang and Stephan ten Brink explaining the differences between OFDM and GFDM.

To sum it up, GFDM uses generalized pulse shaping filters that can span

a subset of $N_\mathrm S$ subcarriers out of a total number $N_\mathrm{S,tot}$ of subcarriers and
a set of $N_\mathrm Y$ symbols

OFDM is a special case of GFDM where $N_\mathrm Y = 1$, $N_S = N_\mathrm{S,tot}$, and the pulse shaping filter is rectangular in frequency domain. This can be implemented using an IFFT at the transmitter and an FFT at the receiver.

The implementation of GFDM is generally not that simple. I could think of an implementation using $N_\mathrm{S,tot}/N_\mathrm S$ IFFTs of length $N_\mathrm{S,tot}$ where at every IFFT just a subset of $N_\mathrm S$ inputs is used and the others are zero padded to realize the frequency allocation of one subset of subcarriers. Each IFFT output would be followd by a specific pulse shaping filter spanning $N_\mathrm Y$ symbols. Finally, all filter outputs would be added up and converted to an analog signal using an digital-to-analog converter (DAC).

The purpose of GFDM is to achieve, amongst others:

lower out-of-band radiation (better spectral efficiency)
higher robustness against carrier frequency offset (CFO)
higher robustness against sampling time offset

when compared to OFDM.

GFDM is also motivated by the fact that in Long-term evolution (LTE) mobile communications, users are allocated so-called physical ressource blocks (PRB). PRBs consist of a set of adjacent subcarriers and symbols.


Comparing two small audio files (e.g., hand clap) Android Dears,
I am new in signal processing and I am dealing with an Android application that should do some audio files (having the same "small" duration) comparison in real time.
My main focus is about comparing two small audio files (about 10 seconds each one) which are related to a hand clap for example (or other sounds related to one short event like ball falling).
In other terms, I should save a reference audio file (or a significant representation you may suggest) and then compare (in real time) a newly recorded audio file with the one already saved.
For instance, if my audio of reference is a hand clap so I decide that the 2 audios are similar if the second one includes a clap even if there are some noises or the 2 sounds intensity is not the same.
I already checked musicg lib for doing this kind of comparison however, I have to convert my recorded audio files to .wav which may add some slowness in my real time application.
I am wondering if there is another alternative that fits better small audio files (about 10 seconds) having the same duration ?
Thanks in advance for your suggestions and help.

What you're interested in sounds like a generic version of a matched filter. I'm not going to be able to speak to the software aspect of how you implement the file handling and "always-on" listening ability of this problem, but from a signal processing perspective, that's what you want to look into. It's worth noting that your filter is unlikely to match any hand clapping or any ball bouncing based on the reference waveform, but the magnitude of the matched filter response should give you a pretty reliable measure of how close the input is. You can use more advanced techniques if you get a lot of false alarms, like requiring the filter response to be a certain SNR by taking constant ambient noise measurements and comparing the amplitude to those.
EDIT: I re-read your question and am convinced this is the correct approach, but I suggest against trying to find a file-level solution to this issue. What you want to do is open the "reference" file and time-reverse the audio data. Then as mic input comes in, you convolve it with the reversed reference waveform (this is matched filtering). A strong response (impulse) from that convolution output indicates a match. In terms of tools you don't even need a powerful audio library, just something that can get the audio in 1-D vector form and then a time-reverse and 1-D convolution (common operations in matrix libraries).


gaussian window for designing filters in DSP I understand that we do not use rectangular window as a window in filter design because of the ringing/gibbs artifacts due to the sidelobes of the resulting sinc function in the time domain.
Consequently, most DSP books mention hanning, hamming windows because they reduce the ringing artifacts.
But how about using a simple gaussian window? It is symmetric, invest DFT will also be gaussian which will have much lesser sidelobes than sinc function, etc.
Is the ringing artifact less for hanning window than the Gaussian window? For me a narrow transition band is not the topmost priority.

First, a Gaussian window will have some parameter controlling the ‘width’ of the Gaussian pulse.  So there is no singular Gaussian window like there would be with many other window types.
Second, because it is a window, it has a defined beginning and end such that the signal function is zero outside the window. The is mathematically equivalent to multiplying your Gaussian function with a rectangular window  this means the resulting frequency response is the convolution of the gaussian frequency response and the rectangular frequency response, which will contain side lobes (ringing). The magnitudes of the side lobes are a function of the window length and the pulse width.  Whether or not the side lobes are better or worse than, say, a Hanning window depends on how you define those parameters and what your metrics are.
Lastly, and this wasn’t explicitly called out but is probably relevant, is discritization. If you are using the resulting window function as part of a discrete time process (in software, for example) the number of samples in the window will effect the frequency response. Every continuous time window function has an infinite bandwidth frequency response because the time domain response is limited. As such, aliasing will occur and the discrete time frequency response will change. The amount it changes is minimized by increasing the number of samples in the window, either by increasing the window length, or increasing the sample rate.


Trying to decompose a signal into sine waves whose wavelengths aren't restricted to 1/N times sample range? I am trying to essentially do a Fourier transform - I want to fit some data with sine/cosine functions. At first I was trying to do this using FFT, but my problem is that the FFT algorithm doesn't seem to provide accurate information about the actual frequencies/wavelengths that are making up the data, i.e., all the wavelengths reported by the FFT are 1/N of the sample window (where N is an integer.)

But, say that the actual wavelength present is some non-integer quotient such as 1/3.5 times the sample window, etc.? Is there a better algorithm than FFT to fit a function with any-wavelength sine waves, instead of only restricting to specific wavelengths?

(I understand that to get a better fit, it probably won't be a fast algorithm like FFT - that is fine if it is a slow method, as long as it can accurately fit to find any arbitrary present frequencies.)

This is my first time posting on this board - let me know if you need any more information about my question. Bonus points if anyone can suggest a way to do this in Python!

There is absolutely no need for explicit interpolation. The DFT can do that for you. Take you signal frame and apply a windowing function to minimise edge artefacts. Afterwards pad the frame with zeros to determine the number to output bins of the DFT. Then apply your FFT and you have an interpolated spectrum with the frequency step size you desire.


Does a filter need to completely attenuate high frequencies to be considered low-pass? I'm looking at a system with a frequency response that is: $$H(e^{j\omega})=\tfrac{1}{4}(2\cos(\omega)+\cos(2\omega)+2)$$

I think the magnitude of this system looks like this (not sure what is the proper name for a function like this):

Would this be considered a low-pass filter even though the output is always greater than zero? If not, what would be the proper name for it?

Ideal low pass filter is supposed to attenuate/remove all the frequencies above the cutoff. but in practice the filters have non zero magnitude even for frequencies in their attenuation region, still they are called low pass.

if the magnitude of filter is quite high for frequencies above cutoff it might be called bad lowpass filter, but is still a low pass because it attenuates the frequencies above cutoff.


Time shifting and time reversal For the problem given in the fig. how to solve it?

My approach:

Let $y(n)$ be the first output (e.g for 1st case output of system $A$ is $y_1(n)$; for 2nd case output of system $B$ is $y_2(n)$, then $$y_1(n)=x_1(-n)$$ $$V_1(n)=y_1(n+1)=x_1(-n+1)$$

And for the 2nd case:

$$y_2(n)=x_2(n+1)$$ $$V_2(n)=y_2(-n)=x_2(-n-1)$$

So if $x_1(n)=x_2(n)=x(n)$, then $$V_1(n)=x(-n+1)$$ $$V_2(n)= x(-n-1)$$

So the answer should be 'C', but the given answer is 'B'.

What is the mistake and how should I proceed?

First of all, since the system A is not LTI, then the ordering of systems is critical in determining the output, and you can weakly conclude that $V_1[n] \neq V_2[n]$ in general, but is it so for all possible inputs? The following is true: if the orders of $N$ LTI systems are interchanged, the resulting outputs would be the same for all possible inputs. But the following is not (necessarily) true: if at least one system were not LTI, then outputs corresponding to different orderings must be different for all possible inputs... It's not true because there can exist some inputs for which the coprresponding outputs could be same in different orderings...

You have solved the two outputs correctly and reached the condition that $V_1[n]=x[-n+1]$ and $V_2[n]=x[-n-1]$. Can you find some inputs $x[n]$  whose transform under these two systems might yield the same output? i.,e., are there any $x[n]$ for which $x[-n-1] = x[-n+1]$ for all $n$ ?

As @loudnoises has exemplified a trivial input of the form $x[n] = 0$ for all $n$ provides a single example to the option b. Indeed all inputs of the form $x[n] = K$ (the constant signal) satisfy the relation that $x[-n+1] = x[-n-1]$ hence there exists infinitely many (for each different K) inputs of the form $x[n]=K$ for which $V_1[n] = V_2[n]$ .


How does transformation matrix and camera matrix change when indexing pixels from -1 to 1 instead of 0 to N-1 I'm trying to warp frames for view synthesis. Specifically, given a frame1 and poses of frame1 and frame2, I'm trying to warp frame1 to frame2's pose/view. For this, I'm indexing pixel locations as $[0,M-1]$ and $[0,N-1]$ for y and x directions respectively. My code is working perfectly (except occluded regions, which lead to holes).
I'm looking at this View Synthesis paper. The code released by authors instead index pixels in the range $[-1,1]$. Directly using transformation and camera matrices is giving invalid transformed coordinates. Can you please explain how I should change the transformation and camera matrices so that it works for pixel indices in the range $[-1,1]$?
More Details:
I've implemented the warping in python using the below equations. Given a pixel location $p_1=[x,y]^T$
$$\hat{p}_1 = [x,y,1]^T$$
$$ P_1 = Z K^{-1} \hat{p}_1$$
$$ \hat{P}_1 = [P_1^T,1]^T $$
$$ \hat{P}_2 = T \hat{P}_1 $$
$$ P_2 = \hat{P}_2[0:3] $$
$$ \hat{p}_2 = K P_2 $$
$$ p_2 = \hat{p}_2[0:2] / \hat{p}_2[2] $$
$p_1$ is the point in frame1 and $p_2$ is the corresponding point in frame2. $P_1$ is 3D world point in the view of frame1 and $P_2$ is the same 3D point in the view of frame2. $T$ is the transformation between views of frame1 and frame2. $Z$ is true depth of point $p_1$ in view of frame1. K is the camera intrinsic matrix. hat ($\hat{p})$ is used for projective coordinates.
I compute the $4 \times 4$ transformation matrix $T$ as follows:
$$T = \begin{bmatrix}
R & t \\
0 & 1 \\
\end{bmatrix}$$
Where $R$ is the rotation matrix and $t$ is the translation.
The $3 \times 3$ camera matrix $K$ is given by
$$K = \begin{bmatrix}
\frac{W/2}{tan(hfov/2)} & 0 & W/2 \\
0 & \frac{H/2}{tan(vfov/2)} & H/2 \\
0 & 0 & 1 \\
\end{bmatrix}$$
Where $H=W=240$ are the height and width of the frames, $hfov=60,\ vfov=45$
In short, what changes should I make to matrices $T$ & $K$?

Before attempting to answer this, let's first see what camera transformation matrix does. Consider the 3d point
$$ P_1 = [X_1, Y_1, Z_1]^T $$
Pre-multiplying by camera matrix gives
$$ K \cdot P_1 =
\begin{bmatrix}
\frac{W/2}{tan(hfov/2)} & 0 & W/2 \\
0 & \frac{H/2}{tan(vfov/2)} & H/2 \\
0 & 0 & 1 \\
\end{bmatrix}
\begin{bmatrix}
X_1 \\
Y_1 \\
Z_1
\end{bmatrix}
$$
$$ =
\begin{bmatrix}
\frac{W/2}{tan(hfov/2)} X_1 + W/2 \\
\frac{H/2}{tan(vfov/2)} Y_1 + H/2 \\
Z_1
\end{bmatrix} $$
As we can see, $W/2$ at $K[1,1]$ scales $x$ value and the $W/2$ at $K[1,3]$ adds an offset. Thus, if $x \in [-1,1]$ then by this scaling and offsetting, we get $x \in [0,W]$.
Conversely, when $x \in [0,W]$, we should use the camera matrix as defined earlier. Hence when $x \in [-1,1]$, we shouldn't add the scaling or offsetting. Thus, only the camera matrix needs to be changed to
$$K = \begin{bmatrix}
\frac{1}{tan(hfov/2)} & 0 & 0 \\
0 & \frac{1}{tan(vfov/2)} & 0 \\
0 & 0 & 1 \\
\end{bmatrix}$$


linear phase and generalised linear phase filters I know that for linear phase filters, all frequency components have equal delay times. That is, there is no distortion of signal due to the time delay of frequencies relative to one another.

However, for filters generalised linear phase, where the phase is of the form $a\omega+b$, and $b$ is not zero. Will the same be true for them? Thanks in advance!

Yes. The time delay of real-coefficient linear-phase N-point FIR filter is (N-1)/2 samples. The time delay of complex-coefficient generalized-linear-phase N-point FIR filter is also (N-1)/2 samples.

You can prove this to yourself.  Design a narrowband linear-phase lowpass FIR filter, and plot its group delay. Then multiply that filter's coefficients by a complex exponential whose frequency is one fourth the sample rate, creating a complex-coefficient filter. Now plot the group delay of the complex filter. Compare the two group delay plots.


Calculation of line of sight system gain I'm trying to calculate the overall gain of the transmitter-receiver system for a line-of-sight wireless transmission.

My question is how to calculate the overall gain of the system. Is it as simple as calculating the gain of the transmitter and receiver separately and then adding them together?

You need to multiply the antenna gains, not add them. Specifically, if the free-space loss (attenuation) is $L_{FS}$, the transmitter antenna has gain $G_T$, and the receiver antenna has gain $G_R$, then the total system loss $L$ is $$ L = \frac{L_{FS}}{G_T G_R}. $$ The system gain $G$ is $$G = \frac{1}{L} = G_{FS}G _T G_R, $$ where $G_{FS}$ is the free-space gain.

Of course, if you're doing the calculation in decibels, then the antenna gains are added: $$ G_{dB} = G_{FS,dB} + G_{T,dB} + G_{R,dB}. $$


Finding periodicity of spikes, matlab In my time series data, I have often spikes with regular periods. Sometimes, I get two or three different periodical spike sequences (two electrical noise sources), some with T1 seconds period, some with T2 second period and some irregular and at random times.

Lets say I have two groups of overlapping periodical spikes in my data, one with a period of 8 seconds and another with 6. I have zeroed the data except for the times where there is spikes. So, my data now looks like a train of spikes. How can I find the periodicity of the two and separate the two groups? Cross-correlation didn't help as many are just single sample spikes. I was thinking of FFT. But I can't see anything in the spectra at those low frequencies, nor when I use FFT on original data. My data has 500 Hz. sampling frequency. I do not want to remove them, just to find the periodicity. Thanks.

My Scipy code that generates data, plots it and performs the fft is listed below.  The code generates two "spike trains" at intervals of 6 and 8 seconds.  The following plot shows the signal in time, and its fft (sorry, the fft title got mangled).

The interesting thing to notice is that the frequencies you are looking for are there.  Said frequencies are

1/8 Hz or omega = 1/8/500 = 0.00025 omega and

1/6 Hz or omega = 1/6/500 = 0.00033

Which is where the first two peaks are in the frequency domain.  Further out are many more repeating peaks/aliases because these are essentially really extreme, intermittent square waves. You could try low-pass-filtering to isolate the fundamental frequencies and eliminate other peaks, but I think that's probably a lost cause at a sample rate of 500 Hz.  Trying to get a 0.2 Hz cutoff filter @ 500 Hz seems a bit much.  Just zooming in on the frequency region of interest might be good enough.  The next step could be to understand exactly how these spike trains are going to cause repetition/aliasing in the frequency domain and look for a certain pattern.

Or try to either LPF and downsample a lot and retry the FFT route,  calculate individual DFT bins using straight-up DFT math or use the Goertzel algorithm to calculate DFT bins for suspected frequency locations.  Maybe try a Goertzel algorithm from 0.0 Hz to 0.2 Hz with 0.01 spacing - just as an example.
Good luck - it doesn't seem like a trivial problem.

from pylab import *
close('all') #close previous plots

fs=500 # generate signal
s1 = hstack(( 1,zeros(fs*6-1) )) # 1/6 Hz
s2 = hstack(( 1,zeros(fs*8-1) )) # 1/8 Hz
s1r = tile(s1,(12,))
s2r = tile(s2,(9,))
s = s1r + s2r
print "shape of s is %s"%(str(s.shape))

sfft=fft(s) # execute fft
sfreqs=fftshift(fftfreq(len(sfft)))
print "shape of sfft is %s"%(str(sfft.shape))

figure(); # generate plots
subplot(211); plot( 1./fs * arange(len(s)) , s ) #, 'o' )
title('signal, s');xlabel('sec');ylabel('mag')
subplot(212); plot( sfreqs , (abs(sfft)) ) #, 'o' )
title('fft(s)');xlabel('$\omega$');ylabel('linear mag')
axis([0.0,0.0005,0.0,25.])
show()

UPDATE: Thinking about it some more, the best solution probably doesn't involve signal processing.  You may want to threshold the samples in to zeros and spikes.  Then take the interval between the spike locations 0 and 1 and look for that interval through the rest of the spikes.  If that interval is consistent with all of your spikes, then you've found that series of spikes and should "subtract out" that series from the spikes.  Then keep iteratively looking at intervals of spikes.  In a noise, imperfect system, a good implementation of this will probably be tricky, but it's probably a better solution to this specific problem.


How harmonics happens when you mix two signals? Doing a search in quora about Why does Amplitude Modulation have sideband frequencies? i found an answer and within it, a user answered with the following statement:
"When you try mixing two signals, there will be something know as harmonics that will occur. This means that frqncy X (+) Frqncy Y will not just form frequency Z but some more frequencies that will be embedded into the system."
Can someone clarify and or exemplify the statement?

Amplitude modulation is a linear operation (excluding the carrier in the classical AM) and does not introduce harmonics due to the modulated signal. Sidebands are not harmonics. They are created as the message signal spectrum is shifted up to carrier frequency.
Frequency mixing of a message signal with a sinusoidal carrier, another definition of AM, does not introduce harmonics.
In a practical setting if there are some nonlinearities in the implementation (a deviation from the ideal case) you may have harmonics. Otherwise No.


What should be IQ sample rate at least? As i read, IQ sample doesn't need nyquist criteria. What is the mathematical representation of this result? Why IQ sample doesn't need nyquist frequency? I know that can be set to nyquist frequency but i want to know why not equal nyquist frequency at least.
Thanks for answers..

As Marcus commented, it certainly does follow the Nyquist criterion. However you just need to keep in mind that I/q sampling is complex so each complex sample can be considered 2 real independent samples and Nyquist will still apply. Also take a look at this answer "Complex sampling" can break Nyquist?


Why are equalizers moving average model (FIR) Let the equalizer coefficients obtained, say $w$ and the input to the equalizer is $y$ that is a noise corrupted received signal, $rx$ then the equalized signal, which is the output of the FIR based equalizer, should be equal to the transmitted signal / source signal. This is the working principle. But, we often perform the operation :

equalizer_output = filter(w,1,rx).

Question 1: Most equalizers are modeled as FIR filter. My Question is why? I  do not understand how this operation gives the input or how is this performing deconvolution?

Question 2:  For example considering a single input single output channel, $x$ is a BPSK modulated source signal and there is multipath fading with 3 taps (not Rayleigh), then the transmitted signal is
$$y(t) = h(1).x(t-1) + h(2).x(t-2) + h(3).x(t-3) + x(t)+ w(t)$$ where $w$ is the Additive White Gaussian measurement / channel noise. The above model is FIR model and not any filter. This is not the equalizer. In order to devise an equalizer inorder to get the input $x$, do we again devise a FIR filter as an equalizer whose input would be $y$? $h$ is the channel impulse response but it is not any model as such.

Question 3: When do we model a channel as an IIR model and then the equalizer would be an IIR or FIR filter?

The reason why almost all linear adaptive equalizers are implemented as FIR filters is that FIR filters are always stable and that there exist relatively simple and effective adaptation algorithms. Note that much work has been done on adaptive IIR filters (e.g., this book by Phillip Regalia), but in practice FIR filters are still the preferred option.

Note that no real-world channel can be ideally equalized by an FIR filter, but that's not as bad as it may seem because an FIR filter can approximate the ideal equalizer reasonably well if its filter length is chosen appropriately and if the adaptation algorithm is tuned to the application. Also, perfect equalization of the channel is not always desired if the channel is not only dispersive (i.e., acts a filter) but if it also introduces noise. In this case a combined criterion is used that takes into account the ISI due to the dispersive channel as well as the degradation due to the noise.


Measuring the "frequency stability" of a signal I'm thinking about ways to automatically identify and analyze lots of musical instrument samples, and would like to find the part of the sample that's the most stable in frequency, without any accidental vibrato or anything.  My first thought was to do a sliding FFT and find the "width" of the peaks (maybe fitting a parabola to many of the neighboring points?), since a constant frequency chunk will produce a narrow peak (green) and one wavering in frequency should produce a wider peak (red):

On a spectrogram the red is a slightly wiggly line, and the green is a more straight line.

Time domain is ok, too.  Some kind of "slide the signal past itself and look for sections that interfere with each other in a constant way instead of phasing in and out"?

But I'm not sure how easy or reliable this method would actually be.  Any other ideas?

Your technique might work, provided you window and so on correctly, but it does not really strike me as the right approach for pitch vibrato for a lot of reasons, mostly that the change in pitch is very small compared to other things which could change the "width" of your spectrum (eg, timbre).

I think you'd be better off actually tracking the fundamental pitch in the time domain using a technique like YIN and seeing how much it varies over time.


What is cross-spectral density- CSD? I asked a question earlier but I didn't get any answer for it. So now I am simplifying it: what are Cross-Spectral Density (CSD) and Power-Spectral Sensity (PSD)? What is their application? How can I get them in MATLAB?

  $$S_{kl}(\omega)=\lim_{T\to\infty}\frac{1}{T}E\{Y_k^*(\omega)Y_l(\omega)\} $$
  $$S_{kk}(\omega)=\lim_{T\to\infty}\frac{1}{T}E\{Y_k^*(\omega)Y_k(\omega)\} $$

  $S_{kl}(\omega)$ is the cross-spectral density (CSD) function between
  general signals $y_k(t)$ and $y_l(t)$, $S_{kk}(\omega)$ is the
  power-spectral density (PSD) of signal $y_k(t)$, $Y_k(\omega)$ is the
  finite Fourier transform of signal $y_k(t)$ at frequency $\omega$,
  $Y_k^*(\omega)$ is the complex conjugate of $Y_k(\omega)$, and
  $E\{\cdot\}$ is the expectation operator.

My earlier question was:
What does 'wavelet power spectrum', 'Auto-power spectrum','cross-power spectrum' means in wavelet application? I was studying about mode shape identification with wavelet method and these terms confused me.

Power-Spectral Density is the distribution of power along the frequency axis.
It is generally used for non-finite energy signals (mostly not limited in time signals), who aren't square-summable.
The signal's PSD is the autocorrelation of the signal's Fourier Transform, as stated by the Wiener–Khinchin theorem.
In Matlab:

N = length(S);
F = fft(S);
F = F(1:N/2+1);
PSD = (1/(2*pi*N)) * abs(F).^2;
PSD(2:end-1) = 2*PSD(2:end-1);
freq = 0:(2*pi)/N:pi;

see: https://de.mathworks.com/help/signal/ug/power-spectral-density-estimates-using-fft.html

Cross-Spectral Density is the same, but using cross-correlation, so you can find the power shared by a given frequency for the two signals using its squared module, and the phase shift between the two signals at that frequency using its argument.

Cross-Spectral Density can be used to identify the frequency response of a noisy LTI system : if the noise is not correlated to the input or output of the system, its frequency response can be found from the CSD of the input and output.


Problem with the definition of Linearity From the high school maths we know that y=mx+c is a linear equation. However, in DSP the linear system must satisfy Additivity properties which y=mx+c does not hold because of +c.
So, is the definition of Linearity different in DSP and Maths? if so why?
Thanks.

Yes, in some sense the definitions are different. I'll give you two points of view, the first one will support your keen observation, the other will provide evidence for the contrary. These two do not conflict with each other, it's a matter of semantics. If the second one confuses you stick with the first one. Prelude concluded, here goes.

Problem Definition 1: $y = mx + c $ means $f(x)=mx+c$

Here we take the usual point of view in which $y$ is assumed to be the output of some manipulations on $x$. We call that something a function, and we can write down the same expression with a little more mathematical elegance: $f(x)=mx+c$. It is now clear that $f(x)$ is in some sense on output of some mathematical manipulation to which $x$ is the input.

Let's refresh criteria for linearity. A function $g(x)$ is linear is it satisfies both of the following conditions:

$g(a+b) = g(a)+g(b)$ for all $a$ and $b$
$g(cx) = cg(x)$ for all constants $c$

Clearly, our favorite function $f(x)$ does not satisfy either of these properties. So yes, from this perspective $f(x)$ is not a linear function. The closest thing to "linear" we can call it is "affine".

Q.E.D.

You may now brace yourself for part 2 of the answer.

Problem Definition 2: $y = mx + c $ means $L(x,y)=y-mx$

Let's take it one step at a time. Suppose you're trying to solve a system of two linear equations. How do you do it? One way is to write down the equations as follows:

$$ \begin{align} y &= m_1x + c_1\\y &= m_2x + c_2 \end{align}$$

Surely, that's how we've all been doing it since seventh grade. Now all you have to do is solve it by substitution or whichever way you prefer. But what do you do when you have a system of equations of more that two variables? Will you write it down like so?

$$ \begin{align} y &= a_1x + c_1z + d_1\\y &= a_2x + c_2z + d_2\\y &= a_3x + c_3z + d_3 \end{align}$$

That doesn't really look right. And for a very good reason. There are many ways of interpreting functions of any number of variables, and it's not just the semantics that's different. To digress for a moment, take the equation $x^2+y^2=r^2$. Almost anyone (visiting this forum that is) will immediately identify it as an equation of a circle. But recall definition of a function!

If we interpret it as $f(x) = \pm \sqrt{r^2-x^2}$ we get two solutions: an upper half of a circle and a lower half of a circle. The entire circle cannot be a solution because it violates the property that in functions, for every input there's at most one unique output.

If we on the other hand interpret it as $f(x,y)=r^2$, we get back the entire circle as a solution, because we're viewing it as a function of two variables equal to a constant. In other words, even though we're written the same expression $x^2+y^2=r^2$, we must define what we're talking about. Otherwise this problem is not well-defined. In one interpretation it is a function $f:\mathbb{R} \rightarrow \mathbb{R}$, in another interpretation it is a function $f:\mathbb{R^2} \rightarrow \mathbb{R}$. Remember all that mumble about domains and ranges in high school? Yeah, this is exactly what it is. Now, back to our myterious topic of linear functions.

Hopefully, by this point you're already has your aha! moment. If not, here's our finish straight. Remember that system of three equations that didn't quite look right? First of all note that it looks affine, because in addition to variables $x$ and $z$ there are constants $d$ as well. Now a nicer way of writing down this system of equations is like so:

$$ \begin{align} -a_1x + y + -c_1z &= d_1 \\ -a_2x + y + -c_2z &= d_2 \\ -a_3x + y + -c_3z &= d_3\end{align}$$

Now we're getting somewhere. As you can see, we can write it out in matrix form as follows:

$$\begin{bmatrix} -a_1 &1 &-c_1 \\-a_2 &1 &-c_2 \\-a_3 &1 &-c_3 \end{bmatrix}\begin{bmatrix} x \\ y \\ z \end{bmatrix} = \begin{bmatrix} d_1 \\ d_2 \\ d_3 \end{bmatrix}$$

Clearly, this is a linear system of equations. Where's the catch? Well, at first it looked like a system of three functions of the form $f:\mathbb{R^2} \rightarrow \mathbb{R}$, and now we're representing it as a single function of the form $f:\mathbb{R^3} \rightarrow \mathbb{R^3}$.

To clarify, this is a single function which takes in a vector in $\mathbb{R}^3$ and returns another vector in $\mathbb{R}^3$. Let's call this function $L(x,y,z)$, precicely $L:\mathbb{R^3} \rightarrow \mathbb{R^3}$. I will let you check that this function is linear. Concretely, if $\begin{bmatrix} a_{11} &a_{12} &a_{13} \\a_{21} &a_{22} &a_{23} \\a_{31} &a_{32} &a_{33} \end{bmatrix}\begin{bmatrix} x \\ y \\ z \end{bmatrix} = \begin{bmatrix} d_1 \\ d_2 \\ d_3 \end{bmatrix}$ and $\begin{bmatrix} a_{11} &a_{12} &a_{13} \\a_{21} &a_{22} &a_{23} \\a_{31} &a_{32} &a_{33} \end{bmatrix}\begin{bmatrix} \alpha \\ \beta \\ \gamma \end{bmatrix} = \begin{bmatrix} \delta_1 \\ \delta_2 \\ \delta_3 \end{bmatrix}$, then

$\begin{bmatrix} a_{11} &a_{12} &a_{13} \\a_{21} &a_{22} &a_{23} \\a_{31} &a_{32} &a_{33} \end{bmatrix}\begin{bmatrix} x+\alpha \\ y+\beta \\ z+\gamma \end{bmatrix} = \begin{bmatrix} d_1 + \delta_1 \\ d_2 + \delta_2 \\ d_3 + \delta_3 \end{bmatrix}$, and
$\begin{bmatrix} a_{11} &a_{12} &a_{13} \\a_{21} &a_{22} &a_{23} \\a_{31} &a_{32} &a_{33} \end{bmatrix}\begin{bmatrix} kx \\ ky \\ kz \end{bmatrix} = \begin{bmatrix} kd_1 \\ kd_2 \\ kd_3 \end{bmatrix}$

In other words (and yes, this is the real reason mathematicians keep constantly coming up with new concise notation!), let $\vec{u},\vec{v} \in \mathbb{R}^3$ ($\vec{u}$ and $\vec{v}$ and 3-dimensional vectors of real numbers). Then

$L(\vec{u}+\vec{v}) = L(\vec{u})+L(\vec{v})$
$L(k\vec{u}) = kL(\vec{u})$

Linear! Q.E.D.

In conclusion, we've explored mysterious subtleties of mathematics of functions and in particular the importance of defining problems well. The function $f(x) = mx + c$ is obviously non-linear (or more precisely affine), and the function $g(x,y) = y - mx$ is linear.

Come back for more interesting stuff. We like giving twisted answers to simple questions.


What should be the correct scaling for PSD calculation using $\tt fft$ I would to calculate the PSD of a signal using FFT however the result do not match with periodogram command.  What  did was as follow :

len=length(signal);
NFFT=2^10;
signal_fft=fftshift(fft(signal,NFFT));
px=signal_fft.*conj(signal_fft)/(NFFT*len);
f=-[NFFT/2:NFFT/2-1]/NFFT*Fs;
plot(f,px);

However I am wondering :

Should I use Fs instead of NFFT in the denominator for calculation signal_fft? Since most websites mentioned it should be Fs, but according to this link the number of FFT points are used. However NFFT makes more sense since that is number of points I have.

The reason is that if I take 1024 point FFT the resolution is different from $2^{20}$ points FFT so NFFT should be used somewhere to adjust for this change.

There is only one correct way of scaling DFT when calculating PSD with RMS values. Given input signal $x$ and its DFT $X$, the exact formula is:

$$\mathrm{PSD}=\frac{2\cdot \hat{X}}{f_s\cdot S} $$

where:

$\hat{X}=|X|^2=X\cdot X^*$ - squared spectrum magnitude

$S=\sum_{i=1}^{N} w_{i}^{2}$ - scaling factor defined as sum of squared samples of window function

$f_s$ - sampling frequency

Scaling coefficient $2$ takes into account removal of energy at negative frequencies (we drop this side of PSD).

Sampling frequency $f_s$ happens to be in the equation since it is a part of the Effective Noise Bandwidth, and PSD is simply PS (Power Spectrum) divided by ENBW. So in the end,  the unit of PSD is $x^2/\mathrm{Hz}$, where $x$ is the quantity that you want to measure, i.e. Volts.

Scaling factor $S$ is compensating for any loss of energy when applying a window function. Obviously for rectangular window (or no window in fact) this doesn't change anything and one, might simply use $N$ instead (length of input sequence).

You did also mention that you were planning to use high number of DFT points ($2^{20}$) since it "makes more sense". Well, it doesn't. If your input sequence has length $N$ (assume 1024), then you can't get better frequency resolution than:

$$\Delta f = \frac{f_s}{1024}$$

Even if you increase number of DFT points to $2^{20}$ by zero-padding, those spectral samples in between will be only interpolated window spectrum - useless garbage. For more detailed explanation please see my other answer.

Another thing to remember. Zero-padding doesn't change the scaling factor! Samples are added after the input sequence had been multiplied by the window, hence the total energy will not change.

Since you wanted to match the calculation with MATLAB's periodogram function, here comes the code. It is a modified version of this one. I tried to comment the code.

% Parameters and signal
Fs = 1000;
dt = 1/Fs;
t = 0:dt:1-dt;
x = cos(2*pi*100*t)';
N = length(x);
% Add a bit of noise
x = x + 0.001*randn(N, 1);
% Use window function of your choice. Rectangular window means no window.
win = rectwin(N);
% The actual number of points used when calculating the DFT might be higher
% this can help in speeding up the computation (next power of 2), but won't
% increase the frequency resolution (only frequency spacing is affected).
% See: https://dsp.stackexchange.com/a/32161/8202
NFFT = N; % N=2^20
xw = x.*win; % Apply windowing here (periodogram does it on it's own)
xdft = fft(xw, NFFT);
psdx = xdft.*conj(xdft); % psdx = abs(xdft).^2
psdx = psdx(1:NFFT/2+1); % Take only first half of the PSD
% Calculate the scaling factor.
% Note that in case of rectangular window this is simply N.
S = sum(win.^2);
psdx = (2/(Fs*S)) * psdx;
freq = 0:Fs/NFFT:Fs/2;
% Calculate the corresponding PSD using MATLAB's in-built
[pxx, ~] = periodogram(x, rectwin(N), NFFT, Fs);

% Plot the PSD's
plot(freq, 10*log10(psdx), 'r')
grid on
hold on
plot(freq, 10*log10(pxx), 'b--')
title('PSD calculation comparison')
xlabel('Frequency [Hz]')
ylabel('Power/Frequency [dB/Hz]')
legend({'DFT', 'periodogram'})

% Plot error between PSD's
figure
diff = pxx - psdx;
plot(freq, diff)
grid on
xlabel('Frequency [Hz]')
ylabel('Difference')
title('Difference between PSDs')

Please note the perfect match between both PSD's:

And here is the difference:

Probably you can notice that there is a bigger error for 1'st and last bin (Nyquist in fact). That is because we multiplied all frequency bins by scaling factor $2$ to compensate for negative frequencies. This shouldn't be done for DC component (1'st bin) and for Nyquist bin (if it's present) - both of them are unique.

But in a real life, Nyquist bin will be removed by anti-aliasing filter and people will remove mean from the signal anyway. Also DC scaling for RMS values is pointless, especially when doing windowing.


Image not divisable by 8? (JPEG Compression) I have learned that an image is to be divided by 8*8 blocks, and the 2D DCT is then applied to these blocks. My question is, what if the image is not divisible by 8? What should I do in this case?

Then it's size is extended by padding zeros or repeating at the end, such that the new size is divisible by 8. Note that reverse repetition at the end reduces the DCT compression size on that block compared to a zero pad approach.

For example if the image was 29 x 27 , then you should add 3 pixels at the right end and 5 lines at the bottom so that it becomes 32 x 32.

Note that you should include the pixel size of the image in the header part so that during decoding process you create the true sized image as 29 x 32


How to index output of convolution product Unfortunately indexing the output of convolution product confuses me. Suppose $x[n]$ is of length $15$ (i.e. $0 <= n <= 14$) and $y[n]$ has the length of $9$ ($0 <= n <= 8$). if $z[n]$ be the output of convolution product of $x[n]$ and $y[n]$, that is $$z[n] = (x * y)[n] = \sum_{k=-\infty}^{+\infty} x[k]y[n-k]$$
As a result, $z[n]$ would has the length of $23$, right?
Now the question is what is the range of $n$ for $z[n]$? Is it $0 <= n <= 22$ or $-4 <= n <= 18$?
From formula for convolution I think the latter is true. But the problem arises when I work with seismic signals.
Suppose $acc[t_i]$ is seismic signal of Kobe earthquake where $0 <= t_i <= 41.99 sec$ with sampling period of $T_s = 0.01 sec$ (length of signal $4200$) as depicted below:

Now it is desired to denoise $acc[t_i]$ using moving average filter with impulse response $h[t_i]$:
$$ h[t_i] =
\begin{cases}
\frac{1}{30},  & \text{if $0 <= t_i <= 0.29$} \\
0, & \text{if $0.3 <= t_i <= 41.99$}
\end{cases}$$
where $T_s = 0.01 sec$. After applying $h[t_i]$ to $acc[t_i]$ by convolution, the output is the sequence $ACC[t_i]$ of length 8399, as shown in following figure.

Obviously the results of filtering process is the interval marked by red bullets. However according to convolution formula, I think I should label time axis between approximately $-21sec$ to $63sec$.
Is that right or something is missing?

If $x[n]$ and $y[n]$ are both causal and starting at index $0$, then the result of convolution will also be causal and it will start at index $0$. Just plug in $n=-4$ in the expression for $z[n]$, you will find that it will be $0$.
$$z[n] = \sum^{\infty}_{k=-\infty}x[k]y[n-k]$$
First non-zero term in above expression is at $n=0$. Because for $k<0$, $x[k]$ will be $0$, you can re-write the above expression as:
$$z[n] = \sum^{\infty}_{k=0}x[k]y[n-k]$$
In the above expression, put $n<0$, you will see that $y[n-k]$ will always be 0.
There will be a transient delay of $N-1$ samples, when you are using Moving-Average FIR Filter of length $N$, not an advance. The delay is due to the fact that you need $N-1$ previous samples of the input before you can produce an average of $N$ samples. You cannot have an output before you feed input to a causal FIR filter. Expecting an output at $-21sec$, is expecting to look into future which is not correct here.
The figures you have posted is showing that you have 8399 samples of output starting at $n=0$. As you can see, approximately at $5sec$, the seismic activity starts, Moving average just smoothens-out the original input signal with a transience of $0.3sec$.


BER performance of a space-time code over Rayleigh channel model For a given space-time code, and in general terms, how would one proceed in the evaluation of the BER performance of this code ?

In order to be more clear, say i have 16 code matrices(who's elements belong to QPSK modulation) of size 2x2 that could be sent over a Rayleigh flat fading channel by a  MIMO system with 2 transmit antennas and 2 receive antennas.
How would one proceed in the evaluation of the BER of such system ?

Thank you!

The classic approach is to evaluate the union bound of the pairwise error probability, taken over the code matrices.

If $X_0$ and $X$ are two different code matrices, then their pairwise error probability is the probability that the receiver mistakes $X$ for $X_0$, when $X_0$ is transmitted. It is written $P(X_0 \rightarrow X)$.

Assuming all code matrices are transmitted with equal probability, and that the code is $\mathcal{X}$, then the average error rate is upper bounded by
$$ P_{\text{err}} \leq \frac{1}{|\mathcal{X}|} \sum_{X_n \in \mathcal{X}} \sum_{X_k \in \mathcal{X},k \neq x} P(X_n \rightarrow X_k) $$
which is just an application of the union bound.

The calculation of BER from $P_{\text{err}}$ for specific channels can get complicated and in fact in many cases no closed form formula can be found. I would suggest the book "Space-Time Block Coding for Wireless Communications" by Larsson and Stoica for more details, in particular section 4, "Error Probability Analysis".


What Makes a Signal Processing Book Good or Less Good? I've recently been disappointed by a few books that have "using MATLAB" in their title, which makes me wonder why. I like MATLAB so it isn't MATLAB. There are other books that I will keep at considerable inconvenience to myself even if I haven't used them for years, and that also makes me wonder why? I'm curious if others have similar opinions and what they might be.

So I'll just name one good book and why.

LINEAR ESTIMATION T.Kailath, A. H. Sayed, and B. Hassibi

The subjects are explained in a clear and insightful way. One feels that they understand when you get to end of the chapter

Then we try the problems and stare at them in blank confusion. Very little is obvious. We now know we don't understand anything.  There is little that a web search will clarify. The problems only yield after thinking and thinking hard.


Estimation of Glottal Closure Instant Glottal Closure Instants (GCIs) are referred to the instances of signiﬁcant excitation of the vocal tract.

I have seen in some papers that I can measure instant of glottal closure in a voiced segment directly from speech waveforms.

Can someone provide me information how detect GCIs ?

How Glottal Closure Instants can be estimated ?

There are many ways to estimate the instant of glottal closure (GCI).  Here (PDF link) is a recent journal article that reviews many of them.  The freely-available VOICEBOX toolbox for Matlab includes a couple GCI estimation techniques.

For most approaches, the fundamental idea is the same: the speech signal is can be viewed as an excitation signal resonating an all-pole filter, and the location of the excitation instants is the instant of glottal closure.

A simple GCI estimation algorithm might work like this.  The all-pole filter (encapsulating the contributions of the vocal tract and glottal source) can be estimated using linear prediction.  Inverse filtering the speech signal with this (time-varying) filter will give the linear predictive residual signal, which is an estimate of the excitation signal.  For normal voices, this signal will look like a noise impulse train.  Peak-picking this signal will give good results for these voices during sustained vowels.  Voice offsets and onsets will be more troublesome, as well as voices that are very breathy or irregular.

(A quick note about your figure: the signal marked as the residual seems to include glottal source signal, and so it not a linear prediction residual in the same sense as what I refer to in the above.)


DDS - Supersampling above the FPGA clock frequency I'm using a Xilinx Virtex 6 FPGA with 200 MHz clock connected to a 1.6 GSa/s DAC to generate sine waves of up to 800 MHz using 8 parallel Xilinx DDS Compiler IP cores.

Initially, I calculate the 32 bit wide phase increment using the formula provided in the manual with a 32-bit wide phase accumulator, 200MHz clock frequency and desired output frequency.

Then 8 DDS modules are instantiated each with a phase offset of $0,\frac 18 \Delta_\phi,\frac28 \Delta_\phi,\ldots, \frac78 \Delta_\phi$. These output 8 samples of a sine wave on every 200 MHz rising clock edge, matching the higher sample rate of the DAC.

The DDS produces correctly sampled output waves up to a frequency of 200MHz. However for f > 200 MHz, the output waves are the wrong frequency. For example, if I input $f =300\,\text{MHz}$, the generated wave is 100 MHz; if I input $f = 400\,\text{MHz}$, the generated wave is 200 MHz.

My suspicions are there is an overflow at some point in the calculation of the phase increment and phase offset values due to their 32-bit width. My question is how would I get around this to correctly produce waves of up to 800 MHz?

I am not concerned with effects of undersampling the higher frequency waves, i just want to see the correct amount of samples per cycle i.e. 2 in the case of $f =800\,\text{MHz}$.

So, as far as I understand your system, you've basically taken your target sample stream and decomposed it into 8 polyphase components; thus, each DDS runs at $\frac18$ of the target 1600 MS/s rate, i.e. at 200 MS/s.

Thus, the phase increment that a single DDS does per sample it produces is 8 times the phase increment of the interleaved samples.

Thus, if you wanted to produce a wave at 200 MHz, that means one period of the output would take 8 samples; therefore the output's phase difference between successive samples is $\frac{2\pi}8=\frac\pi4$. And thus, the period, from the point of view of the individiual DDS has a period of 1 sample, meaning $\Delta_\phi =0$! That makes sense: if you want to produce an exact 200 MHz wave, your 8 polyphase components would always be the same samples of a sine wave period.

So, your behaviour is fully expected. When your output frequency is 200 MHz, the individual DDS'es frequency must be 0. You basically work based on aliases.

This has nothing to do with integer overflows, but is just like the math for polyphase signal generation is.

So, nothing wrong here, you'll just have to figure out what the right phase offset per DDS is.


Matched filter loss due to frequency offset I'm working with a pulsed radar signal which uses linear chirp pulses. I perform pulse compression by applying a matched filter of the simulated chirp waveform and recording the location and power of the peak. Often times the recorded signal has a Doppler shift, so there is some inherent loss in the pulse compression process. This value is not large, but I need to quantify this expected loss so I can properly measure the actual received power of my signal.

How can I calculate the expected compression loss at the peak due to a frequency offset?

PS, I understand that the location of the peak will also move due to range-doppler coupling, but I understand that phenomenon already.

Short answer; I'll try to come back and improve this later.

Recall that a matched filter is equivalent to a sliding correlator. At each time instant, the sliding correlator multiplies the last $T$ (where $T$ is the pulse duration) of the input signal by the conjugate of the pulse shape and coherently integrates the result over the pulse duration. If there is no frequency offset, then the result after multiplying by the conjugated pulse shape is a constant (DC) value, so you don't observe any loss.

If there is a frequency error, then the result after multiplying by the conjugate of the pulse shape is a sinusoid whose frequency is commensurate with the frequency offset. The scalloping loss is incurred is due to the fact that you're integrating the value of said sinusoid over a time window of length $T$. Therefore, the amount of attenuation can be well-modeled by the frequency response of a boxcar filter of length $T$:

$$
\text{rect}\left(\frac{t}{T}\right) \Leftrightarrow T \text{ sinc}\left(\frac{2\pi f T}{2}\right)
$$

You can make the following observations from this relationship:

As the pulse duration $T$ increases, the frequency response of the correlator rolls off more quickly; for a given frequency offset $f$, you'll experience more attenuation from the ideal peak value.
As the frequency offset $f$ increases, the frequency response of the correlator rolls off more quickly as well.

These combine to tell you what you probably already intuitively know; the longer you coherently integrate a signal (i.e. the longer the pulse duration in your case), the more sensitive the result will be to frequency offset.


Why Zero Padding in the Center of the DFT Interpolates / Upsamples the Signal (Sinc Interpolation / DFT Interpolation / Periodic Interpolation) I'm experimenting with the Inverse Discrete Fourier Transform. Starting from the two-cycles continuous $x(t)$ signal below:

I have the discrete signal $x(n) = \{ 1, 0, -1, 0, 1, 0, -1, 0 \}$ leading to the 8 points DFT $X_0(n) = \{ 0, 0, 4, 0, 0, 0, 4, 0 \}$

Now, if I use the IDFT on $X_0$, I obtain $x_0$ looking like that (the blue curve is the real part of the IDFT $Re[x_0(t)]$):

Here $x_0(n) = x(n)$ for integer values ($n = 0, 1, 2, \dots 7$). I understand I do not get back the continuous $x(t)$ function because of aliasing. I would explain that by saying on $x_0(t)$, there is a second signal, above the Nyquist frequency, and "riding" the "carrier"1

I read about zero padding, so I tried to add $k$ extra zeros in the middle of $X_0(t)$ and now, if I perform the IDFT, I obtain the following results:

1. With 8 extra values

I have $X_{k=8}(t) = \{ 0, 0, 4, 0, \underbrace{\mathbf{0, 0, 0, 0, 0, 0, 0, 0}}_\text{8 extra bins},  0, 0, 4, 0 \}$, leading to $x_{k=8}(t)$:

2. With 24 extra values

I have $X_{k=24}(t) = \{ 0, 0, 4, 0, \underbrace{\mathbf{0, 0, 0, 0, \dots, 0, 0, 0}}_\text{24 extra bins},  0, 0, 4, 0 \}$, leading to $x_{k=24}(t)$:

3. The problem

Adding more bins proportionally decreased the amplitude of the rebuild signal and increase the frequency of the signal riding the carrier. I think I understand both phenomenons.

However, I can't find an intuitive way of explaining why, at integer positions (the orange dots), $x_k(n)$ reproduces more and more accurately the original shape of the $x(t)$ signal.

1Do not hesitate to edit the question if I don't use the correct vocabulary here.

What you are experiencing is technically called interpolation by DFT; i.e., interpolating the time-domain sequence $x[n]$ by properly zero staffing the middle portion of it's DFT $X[k]$ (and taking the inverse DFT to get the time domain interpolated sequence).

Normally the signal interpolation is described and implemented in time-domain. However it's also possible to implement it using its frqeuency domain expression...

The following may be helpful:
Interpolation explained in time-domain:

$$ x[n] \rightarrow ({\uparrow L}) \rightarrow w[n] \rightarrow \boxed{LPF} \rightarrow y[n]$$

Input sequence $x[n]$ is of length $N$, expanded sequence $w[n]$ is of length $M = N \times L$, and the LPF has a Gain = L and discrete-time cutoff frequency of $\omega_c = \pi/L$ radians per sample. the interpolated output sequence $y[n]$ is also of length $M = L \times N$.

Let's look at the relation in the frequency domain :
if $X[k]$ is the N-point DFT of $x[n]$, then $W[k] = X[k]$ is the $M  = L \times N$ point DFT of the sequence $w[n]$ (Note: $W[k]$ is indeed an L-fold copy of $X[k]$.). Let the DFT of LPF be $H[k]$ which is also $M = L \times N$ point.

$$ H[k] = \begin{cases}{ L ~~~,~~~-N/2 \leq k \leq N/2 \\ 0 ~~~~, ~~~~ \text{otherwise} }\end{cases} $$

Then the DFT of the interpolation output $y[n]$ is
$$Y[k] = H[k] W[k] $$

After the multiplication $Y[k]$ becomes :
$$ Y[k] = \begin{cases}{ L X[k] ~~~,~~~-N/2 \leq k \leq N/2 \\ 0 ~~~~, ~~~~ \text{otherwise} }\end{cases} $$

This is exactly what you are implementing in your zero staffed DFT of $X[k]$. You try to obtain $Y]k]$ by staffing the middle portion of $X[k]$ to make it length $M$. And you can also see why your magnitude is missing by $1/L$ as you did not multiply it by $L$.

The implicitly  implemented lowpass filter refers to how you obtain $Y[k]$, by zero staffing $X[k]$ into length $M$, instead of obtaining it by the explicit product $Y[k]= H[k] X[k]$ which shows the filtering in freqency domain explicitly.

The reduction in magnitude can be explained both as due to the increase in the length of the interpolated sequence (hence inverse DFT scaling) or as due to the (missing) gain of the implicit interpolation lowpass filter that is implied by zero staffing the DFT. To correct the amplitude mismatch, simply multiply the interpolated sequence by the interpolation factor.


Haar series approximation of a function? I'm trying to code the Haar wavelet approximation of a function.
I'm a complete newbie.

Edit
Below is a pseudo-code of my algorithm for $f: [0,10] \to \mathbb R$, $f(x)=x$.
I picked it here and there on google.

1) Haar wave formula of order $N$  for $x \in [0, \infty)$
$$
\psi_{N,j}(x) = \frac{1}{\sqrt{N}} \cdot\begin{cases}
1 & j =0 \\
+2^{p/2} & a \le x - \left \lfloor \frac{x-a}{L}\right\rfloor L< b \\
-2^{p/2} & b  \le x - \left \lfloor \frac{x-a}{L} \right\rfloor L < c\\
0 & \text{otherwise}
\end{cases}
$$
where
\begin{align}
a &= \frac{q-1}{2^p},
b = \frac{q-0.5}{2^p},
c = \frac{q}{2^p},
L = c-a, \\
p &= \left\lfloor \log_2(j+1) \right\rfloor,
q = j +1 - p.
\end{align}

2) Haar series coefficients
$$
\hat \psi_{N,j} = \int_0^{10} f(t) \psi_{N,j}(t) dt
$$

3) Haar series
$$
f(x) \approx \sum_{j=0}^N \hat \psi_{N,j} \cdot \psi_{N,j}(x)
$$

4) For $N=4$, this gives the graph

Thanks to Jan's code. I made my implementation work.
The code below compare: Haar vs Fourier vs Chebyshev.

from __future__ import division
from mpmath import *

# --------- Haar wavelet approximation of a function
# algorithm from : http://fourier.eng.hmc.edu/e161/lectures/wavelets/node5.html
# implementation only handle [0,1] for the moment: scaling and wavelet fcts need to be periodice

phi = lambda x : (0 <= x < 1) #scaling fct
psi = lambda x : (0 <= x < .5) - (.5 <= x < 1) #wavelet fct
phi_j_k = lambda x, j, k : 2**(j/2) * phi(2**j * x - k)
psi_j_k = lambda x, j, k : 2**(j/2) * psi(2**j * x - k)

def haar(f, interval, level):
    c0 = quadgl(  lambda t : f(t) * phi_j_k(t, 0, 0), interval  )

    coef = []
    for j in xrange(0, level):
        for k in xrange(0, 2**j):
                djk = quadgl(  lambda t: f(t) * psi_j_k(t, j, k), interval  )
                coef.append( (j, k, djk) )

    return c0, coef

def haarval(haar_coef, x):
    c0, coef = haar_coef
    s = c0 * phi_j_k(x, 0, 0)
    for j, k ,djk in coef:
            s += djk * psi_j_k(x, j, k)
    return s

# --------- to plot an Haar wave
interval = [0, 1]
plot([lambda x : phi_j_k(x,1,1)],interval)

# ---------- main
# below is code to compate : Haar vs Fourier vs Chebyshev

nb_coeff = 5
interval = [0, 1] # haar only handle [0,1] for the moment: scaling and wavelet fcts need to be periodice

fct = lambda x : x

haar_coef = haar(fct, interval, nb_coeff)
haar_series_apx = lambda x : haarval(haar_coef, x)

fourier_coef = fourier(fct, interval, nb_coeff)
fourier_series_apx = lambda x: fourierval(fourier_coef, interval, x)
chebyshev_coef = chebyfit(fct, interval, nb_coeff)
chebyshev_series_apx = lambda x : polyval(chebyshev_coef, x)

print 'fourier %d chebyshev %d haar %d' % ( len(fourier_coef[0]) + len(fourier_coef[1]),len(chebyshev_coef), 1 + len(haar_coef[1]))
print 'error:'
print 'fourier', quadgl(  lambda x : abs( fct(x) - fourier_series_apx(x) ), interval  )
print 'chebyshev', quadgl(  lambda x : abs( fct(x) - chebyshev_series_apx(x) ), interval  )
print 'haar', quadgl(  lambda x : abs( fct(x) - haar_series_apx(x) ), interval  )

plot([fct, fourier_series_apx, chebyshev_series_apx, haar_series_apx], interval)


Signal segmentation into states ( Stop/mov) My question concerns signal segmentation.
My signal is the green one, which represents the distance from the first point to the rest of points.
I would like to cut the  trajectory into segments that each represents  the current status of the person, is he currently making a move, or making a stop?
This is my plot:

the green line is the distance measures smoothed using savitzky golay.
The second line in red represents the results I would like to have, when its =1 its a move, and when it is equal to 0 it is a stop in reality.As you can see, it is somehow easy to segment the trajectory using the trends, when it is increasing or decreasing, we should have a move. When it is constant, we should have a stop.

My question is, what are the methods that I can use and that can segment a signal using the trend change ?
Thank you

It seems like you already got all the right components. Your goal is to detect segments where speed!=0 in the data. Since speed is the derivative of position, you could transform the data into the first difference series - simply subtract the preceding point from each data point, i.e. $x_t - x_{t-1}$, and you'll have a crude estimate of the derivative.
The series of differences then can be segmented with any standard method for peak or changepoint detection. In particular, since you know the background value (zero), and the speed appears roughly constant during each movement episode, it seems that epidemic detection of changes in mean should be able to solve this nicely. I am aware of a couple packages for this in R, but it's also relatively easy to implement it yourself if it's not available for your programming language (just search for "epidemic changepoint detection"). Multivariate detectors in the same framework are also available, if you were interested in incorporating several sensors to detect movement, for example.


How to perform an FFT on a signal with a sampling rate of 44100? This is probably a very simple Signal Processing question but I can't quite work it out. I'm currently writing an android app which needs to perform an fft on audio data taken from the phone's microphone. I'm currently reading in the audio data into a byte array buffer at a sampling frequency of 44100Hz.

 int bufferSize = AudioRecord.getMinBufferSize(44100,
                        AudioFormat.CHANNEL_IN_STEREO, AudioFormat.ENCODING_PCM_16BIT);

                final AudioRecord recorder = new AudioRecord(MediaRecorder.AudioSource.MIC, 44100,
                        AudioFormat.CHANNEL_IN_STEREO, AudioFormat.ENCODING_PCM_16BIT,
                        bufferSize);

                final AudioTrack audioTrack = new AudioTrack(
                        AudioManager.STREAM_MUSIC, 44100,
                        AudioFormat.CHANNEL_IN_STEREO, AudioFormat.ENCODING_PCM_16BIT,
                        bufferSize, AudioTrack.MODE_STREAM);

                final byte[] buffer = new byte[bufferSize];
                final int finalBufferSize = bufferSize;

Once the buffer has been filled I convert the data to an array of Complex numbers and pass the data into an fft function. My problem is that the fft function I have cannot work with 44100 samples as it requires the length of the array of samples to be a power of 2. I've found what I think are some solutions online but I'm a little unsure if I'm fully understanding them correctly.

As I understand it I can either pad the array out with zeroes in order to make the array length a power of 2 or I can take a fixed fft size, e.g. 1024 samples. If I use the second solution do I need to divide the original array of samples into 1024 chunks and run an fft on each? Any advice on how I should go about this would be appreciated.

Here is the fft class I'm using:
http://introcs.cs.princeton.edu/java/97data/FFT.java.html

A 44100 Hz sample rate just means that you get 1 sample approximately every 22.67 microseconds. It also means you get N samples in N*22.67 microseconds, and you're free to choose N. If you choose N=44100, you get 44100 samples in one second. If you choose N=1024, you get 1024 samples in 23.22 milliseconds.


Hilbert transform to compute signal envelope? I've heard that the Hilbert transform can be used to calculate the envelope of a signal.  How does this work? And how is this "Hilbert envelope" different from the envelope one gets by simply rectifying a signal?

I'm interested specifically in finding a way to calculate an envelope for use in dynamic range compression (i.e., "turning down the volume" of the loud parts of an audio signal automatically).

The Hilbert transform is used to calculate the "analytic" signal. See for example http://en.wikipedia.org/wiki/Analytic_signal. If your signal is a sine wave or an modulated sine wave, the magnitude of the analytic signal will indeed look like the envelope. However, the computation of the Hilbert transform is not trivial. Technically it requires a non-causal FIR filter of considerable length so it will require a fair amount of MIPS, memory and latency.

For a broad band signal, it really depends on how you define "envelope" for your specific application. For your application of dynamic range compression you want a metric that is well correlated with the the perception of loudness over time. The Hilbert Transform is not the right tool for that.

A better option would be to apply an A-weighted filter (http://en.wikipedia.org/wiki/A-weighting) and then do a lossy peak or lossy RMS detector. This will correlate fairly well with perceived loudness over time and is relatively cheap to do.


Time-differentiation property of Fourier transform for $\sin(\omega_0t)$ and $\cos(\omega_0t)$ As I know, the Fourier transform has the below property which is called time-differentiation:

$$
\frac{dx(t)}{dt}\leftrightarrow j\omega X(j\omega)
$$

and the Fourier transform of the cosine and the sine signals are $\pi(\delta(\omega-\omega_0)+ \delta(\omega+\omega_0))$ and $\frac {\pi} {j}(\delta(\omega-\omega_0)- \delta(\omega+\omega_0))$, respectively.

The derivative of $\sin(t)$ is $\cos(t)$. However, the Fourier transform of $\cos(t)$ cannot be obtained from multiplying the Fourier transform of $\sin(t)$ with $j\omega$.

could any one help me.

It does work. Note that the Fourier transforms you stated in your question are not the transforms of $\cos(t)$ and $\sin(t)$, but of $\cos(\omega_0t)$ and $\sin(\omega_0t)$.

Let's Assume we know the transform of $\sin(\omega_0t)$:

$$\mathcal{F}\{\sin(\omega_0t)\}=\frac{\pi}{j}[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)]\tag{1}$$

Now we take the derivative of $\sin(\omega_0t)$:

$$[\sin(\omega_0t)]'=\omega_0\cos(\omega_0t)\tag{2}$$

Using $(1)$ and the differentiation property of the Fourier transform, the Fourier transform of $(2)$ must be

$$\mathcal{F}\{\omega_0\cos(\omega_0t)\}=j\omega\cdot \frac{\pi}{j}[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)]\tag{3}$$

Noting that for any function $f(\omega)$ that is continuous at $\omega=\omega_0$ we have

$$f(\omega)\delta(\omega-\omega_0)=f(\omega_0)\delta(\omega-\omega_0)\tag{4}$$

and, consequently, we can rewrite $(3)$ as

$$\begin{align}\mathcal{F}\{\omega_0\cos(\omega_0t)\}&= \pi[\omega_0\delta(\omega-\omega_0)-(-\omega_0)\delta(\omega+\omega_0)]\\&=\omega_0\pi[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)]\end{align}\tag{5}$$

as expected.


amplitude of upsampled and downsampled signal without filter Given:
$$
DTFT\{x[n]\}=X(\omega)= \begin{cases}
1 & |\omega| \leq 2/\pi
\\
0 & 2/ \pi < |\omega| < \pi
\end{cases}\ \ \ \ \ (periodic\ 2\pi)
$$

If I downsample $X(\omega)$ by M. I get:

$$
\begin{aligned}
x_d[n]      &= x[nM] \\
\\
X_d(\omega) &= 1/M\ X(\omega/M)
\end{aligned}
$$

Now what happens if I up-sample $x[n]$ again by $L$?

why can't I just let $M = 1/L$ and substitute into down-sampling equation to get up-sampling equation?

$$
\begin{aligned}
x_i[n]      &= x[n/L] \\
\\
X_i(\omega) &= L\ X(L \omega)
\end{aligned}
$$

I'm a little bit confused because my book is telling me that upsampling equation doesn't scale the amplitude by L:

$$
X_i(\omega) = X(L\omega)
$$

Just wondering what happened to the L term when upsampling?

(page 112, Schaum's Outline, Digital Signal Processing, Second Edition)

If I look up DTFT time scaling property, I find:

$$
DTFT\{\ x[an]\ \} = 1/a X(\omega / a)
$$

Which seems to confirm that upsampler should have an L term ..

Close, but not the right formula for the up-sampler:

$$x_i[n] = x[n/L]$$

The "up sampler" is really an expander that inserts zeros between the samples of x[n], and it has the formula:

$$
x_u[n] =\begin{cases}
x[n/L] & n=0, \pm L, \pm 2L, ... \\
\\
0 & otherwise
\end{cases}
$$

The DTFT of the up-sampler "expander" can be derived as follows:

using $\delta(n)$ we can rewrite $x_u[n]$ as:

$$
x_u[n] = \sum_{k=-\infty}^{\infty}x[k]\delta[n-kL]
$$

Applying definition of DTFT to $x_u[n]$:

$$
\begin{aligned}
X_u(\omega) &= \sum_{n=-\infty}^{\infty} \left(x_u[n] e^{-j\omega n} \right) \\
\\
X_u(\omega) &= \sum_{n=-\infty}^{\infty} \left( \sum_{k=-\infty}^{\infty}x[k]\delta[n-kL] \right) e^{-j\omega n} \\
\end{aligned}
$$

Applying sifting property of $\delta$:

$$
\begin{aligned}
X_u(\omega) &= \sum_{n=-\infty}^{\infty} x[k] e^{-j\omega L k} \\
\\
X_u(\omega) &= X(\omega L)
\end{aligned}
$$

Thus, the gain of up-sampler "expander" in freq domain is unity.


Why does this not work: Alt. way to compute the period of a discrete periodic signal that is a sum of 2 complex exponentials I want to compute the fundamental period of the following discrete time signal:

$ x[n]= \exp^{(j\frac{2\pi}{3})n} + \exp^{(j\frac{3\pi}{4})n}$

I know I can do this by taking the fundamental periods of the first term $\exp^{(j\frac{2\pi}{3})n}$, $N_1=3$ and of the second term $\exp^{(j\frac{3\pi}{4})n}$ $N_2=8$ and taking the least commmon multiple which gives $N=24$.

However I tried to express $x[n]$ as a single term as:

$ x[n]=\exp^{(j\frac{8.5\pi}{12})n} (\exp^{(j\frac{0.5\pi}{12})n}-\exp^{(j\frac{-0.5\pi}{12})n}) $

$ x[n]=2\exp^{(j\frac{8.5\pi}{12})n}.\cos(\frac{0.5 \pi n}{12}) $

Since the modulus of the complex exponential represents the magnitude of the signal, the fundamental period then equates to $N=48$ using only the cosine term. I tried this for varying frequencies and using this method I always got a the value for the period that was twice the correct value.

I know there is a flaw in my reasoning but could someone point that out ? And is there a reason why the fundamental period computed in this way always comes out to be twice that of the correct value ?

The flaw in your reasoning is that you only consider the cosine term. The cosine term has a period of $48$. But there's also the complex exponential, which is not just a magnitude but a function of $n$. You can show that that exponential term is periodic with a period of $48$. However, you cannot conclude that the multiplication of the two terms must also have a period of $48$.

Multiplying the two terms results in a period of $24$ because

$$e^{j\frac{8.5\pi}{12}n}=-e^{j\frac{8.5\pi}{12}(n+24)}$$

and

$$\cos\left(\frac{0.5\pi n}{12}\right)=-\cos\left(\frac{0.5\pi (n+24)}{12}\right)$$

so the sign changes cancel in the product.


Promote the Orthogonality between Rows of $ S $ I have a question.

Suppose we want to solve an optimization problem:

Consider $S \in \mathbb{R}^{N \times T}, T >> S$

$$\min_{S} f(S) \mbox{ s.t. } SS^T \mbox{is diagonal}$$

Which means each rows of the matrix $S$ is mutually orthogonal.

I am suggested to solve this alternative problem by the following method:

$$\min_{S} f(S) + \|\mathcal{P}(SS^T)\|_1$$, in which $\mathcal{P}$ is a projection onto the off-diagonal indexes.

But I don't think this $\ell_1$ penalty will promote the orthogonality between each rows by simply promote the sparsity of the off diagonal elements, since we are not doing any actions like block coordinate descent to promote the orthogonality.

Any suggestions?

First, a warm welcome to StackExchange! Consider this as a bit of brainstorming.

Given that the diagonal elements are positive, I think there is a particular geometry of your parameters: They live in the Generalized Stiefel Manifold, where each point is a scaled orthonormal frame of dimension $T$ embedded in $\mathbb{R}^N$.

Let me rewrite :
$$
\begin{align*}
SS^\top &= \Lambda\\
\frac{1}{\sqrt{\Lambda}}SS^\top\frac{1}{\sqrt{\Lambda}^\top} &= \frac{1}{\sqrt{\Lambda}}\Lambda\frac{1}{\sqrt{\Lambda}^\top}\\
B SS^\top B^\top &= B\Lambda B^\top\\
OO^\top &= I
\end{align*}
$$
where $\Lambda$ is the $N\times N$ diagonal matrix, $B=\frac{1}{\sqrt{\Lambda}}\succ0$ and $O=BS$. Similarly, $S=\sqrt{\Lambda}O$. Division and square roots are element-wise. Now $O$ lives on the usual Stiefel manifold. When $B$ is known, the problem can simply be solved by a Riemannian optimization on the Generalized Stiefel Manifold (see manopt). However, now, this is not the case.

So let's take a second look at the problem at hand:
$$\min_{S} f(S) \mbox{ s.t. } SS^T \mbox{is diagonal}$$

This can be formulated in terms of the new variables:
$$\min_{\Lambda,O} f(\sqrt{\Lambda}O) \mbox{ s.t. } OO^T =I \,\wedge\, \Lambda>0$$

where $O$ is Stiefel, and $\Lambda>0$. The latter can also be denoted by a vector $\mathbf{\lambda}\equiv \text{diag}(\Lambda)$ and $\lambda>0$. The true exponential and logarithm maps for the Stiefel manifold are available (for instance in Manopt) and $\lambda$ has a simpler constraint (for projection). So from this point on, I am being more heuristic and would suggest to alternatively optimize these two:

Fix $\lambda$, optimize $O$ on the Stiefel manifold using the chain rule:
$$
\nabla_O f(\sqrt{\Lambda}O) =  \sqrt{\Lambda} \nabla_{S} f(S)
$$
Fix $O$, optimize for $\sqrt{\lambda}$ (can also use positivity constraint):
$$
\nabla_{\sqrt{\Lambda}} f(\sqrt{\Lambda}O) =  O \nabla_{S} f(S)
$$

You can of course use gradient descent but also other sophisticate algorithms like Riemannian-LBFGS / trust-region.


Determine if $ y[n] = ny[n-1] + x[n]$ is linear time invariant and BIBO stable Check if the following system is linear time invariant and BIBO stable..
$$
y[n] = ny[n-1] + x[n]
$$

for $n\ge 0$. We are also given that the system is at rest (i.e. $y[−1] = 0$).

I know that to show linearity, we show that
$
H[ax_1 + bx_2] = aH[x_1] + bH[x_2]
$
except that the question has a $y[n-1]$ on the RHS which is throwing me off.

Edit: My solution:
$
LHS = $ $$H[ax_1[n] + bx_2[n]]
$$$$
= n(ay_1[n−1] + by_2[n−1]) + (ax_1[n]+bx_2[n])
$$$
 RHS =$$$  aH[x_1] + bH[x_2]
$$$$
=  any_1[n−1]+ax_1[n]+bny_2[n−1]+bx_2[n]
$$
So $
 RHS = LHS
$

Edit 2: New issue when I try and show time invariance.

If $ x_2[n] = x1[n-k] $ then prove $y_2[n]  = y_1[n-k]$.

BUT $RHS= $
$$
(n-k)y_1[n-k-1] + x_1[n-k]
$$ And $LHS = $
$$ny_2[n-1] + x_2[n] = ny_2[n-1] + x_1[n-k]$$

I'll give you some hints that hopefully will allow you to do your homework yourself.

In order to check linearity, write down $3$ difference equations: one for the response $y_1[n]$ to an input signal $x_1[n]$, one for the response $y_2[n]$ to an input signal $x_2[n]$, and one for the response $y_3[n]$ to an input signal $ax_1[n]+bx_2[n]$. Multiply the first difference equation by $a$, the second one by $b$, and add them. If the resulting difference equation for the sequence $ay_1[n]+by_2[n]$ is the same as for $y_3[n]$, then you can conclude $y_3[n]=ay_1[n]+by_2[n]$, and, consequently, the system must be linear.

For checking time-invariance, write down the difference equation for the response $y_2[n]$ to the input $x[n-k]$. Then replace in the original difference equation the index $n$ by $n-k$ and check if the two difference equations are the same. If they are, the system is time-invariant.

To answer the question about BIBO stability, ask yourself what the output signal will be for $x[n]=\delta[n]$ and draw your conclusion.

Now that you've figured out the solution I would like to point out that you can also find an explicit (non-recursive) expression for $y[n]$:

$$\begin{align} y[0]&=x[0]\\
y[1]&=1\cdot x[0]+x[1]\\
y[2]&=2\cdot 1\cdot x[0]+2\cdot x[1]+x[2]\\
y[3]&=3\cdot 2\cdot 1\cdot x[0]+3\cdot 2\cdot x[1]+3\cdot x[2]+x[3]\\
&\vdots\\
y[n]&=\sum_{k=0}^n\frac{n!}{k!}x[k],\qquad n\ge 0\tag{1}
\end{align}$$

from which it is very straightforward to show that the system is linear and time-varying. Furthermore, from $(1)$ it should be clear that the system is not BIBO-stable.


Measuring Phase Difference at DC bin aliased frequency using FFT I have a code in which two complex signals (I/Q) are created having same frequency but with phase difference. FFT of both signals is taken and atan2 is used to compute phase after which the difference between both phases is calculated.

While changing the input frequency it was observed that for the case where the signal frequency is integer multiple of Fs such that its FFT peak lies on DC bin, the correct phase difference is still obtained. How can phase calculated from DC bins give the correct value. Can anyone explain the underlying logic?

The "DC bin" does not represent DC exclusively, but has also some low frequency energy. How much depends on the fft size and sampling rate. Simple example: if your fft size is 1000 and your sampling rate is 1000Hz, then the dc bin represents frequencies from 0 to 0.5Hz, as the distance between fft bins is 1Hz. So, a phase can be calculated for the DC bin.


Serial filtering using FFT lets say I have 3 audio filters for a realtime DSP application. For simplictity each has length 256 (as well as the input signal). The filters should work in series. Starting with filter IR h1(n) and input x(n), zeropadding to 512, FFT, then multiplication and IFFT, delivers the first stage output Y1(n).

Since they input vector was padded to 512 the output has 512, too.

Second filter stage: input to this stage is y1(n) which has already 512 samples, and has to be convolved with filter h2(L=256). So y1 and h2 have to be zeropadded to 788, FFT, Multiplication, IFFT... Output y2 has lenght 768

Third stage: will add another padding, so final output vector will have length 1024

Since my app is realtime I would, for performance reason, like to avoid multiple FFT->iFFT for each stage. Moreover there will be a growing input length for the FFT with each filter stage.
I would prefer doing FFT for each filter individually to get H1,H2,H3 then do X*H1*H2*H3, then IFFT. As far as I understood from a previous post, this would work, if I would zero pad each filter (and input) already in the beginning to length L=N+N1+N2+N3. Unfortunately this I want to avoid since some filters are laying already pre-FFt´ed in big buffers already from the app start to grant fast access, and others (Equalizers e.g.) can be optionally switched on/off during runtime. So I have no way to know at program start to know how much "padded" zeros I will need to apply to those pre-FFT´ed filters.

Is there a way to get arround this? Maybe zero padding those static filters in frequency domain "on the fly" if there should be added e.g. another EQ filter on runtime? If yes, how?

No, if you do zero-padding of $H_1$ in frequency domain.

But yes if you do zero-padding of $h_1$ in time domain.
Basically you should have two FFT versions of $H_1$ ($H_{1,512}$ and $H_{1,1024}$). Because $h_1$ is static, you can precompute both, or you can compute $H_{1,1024}$ only when required.

Let $a = IFFT(FFT(x) \times H_{1,512})$ and $b = IFFT(FFT(x) \times H_{1,1024} \times H_2
 \times H_3)$. Size of a is 512 and size of b is 1024. For example, using overlap add method, we can freely change fft-size between 512 and 1024, but it still process 256-input per fft.

===a====
     +
    ===a====
         +
        ===b============
             +
            ===b============
                 +
                ===a====


Estimation Theory - Basic Question I'm just starting estimation theory for my DSP course and I'm a bit confused about the classic example given everywhere which aims to show that two unbiased estimators can have different variances.

I'm just wondering how they did that step that I have a question mark beside. I understand intuiviely that the average of x[n] will be A but I can't mathematically figure that out with the E[] and also what happens to 1/N. If someone could just explain that step in detail.

Another Question: What's the difference between A(hat) and E[A(hat)], I mean shouldn't A(hat) give you A anyways? Then what's the point of the E[A(hat)]?

The expectation of $x[n]$ is

$$E(x[n])=E(A+w[n])=E(A)+E(w[n])=A+0=A\tag{1}$$

because $A$ is deterministic and the noise $w[n]$ is assumed to have a mean of zero. So in that formula you sum up $N$ terms that are equal to $A$, and then you divide that sum by $N$, leaving you with the value $A$.

As for the expectation of the estimated value $\hat{A}$, note that $\hat{A}$ is a random variable because it necessarily depends on the data $x[n]$, which are random, so computing its expected value (and its variance) makes sense.


How to reduce frequency resolution for high sampling rate and lot of samples? I'm kinda new in this field, and I have a following situation. There is a 1s long audio signal, with Fs = 96000Hz, and I have 96000 samples accordingly. If I do a single FFT on the entire signal, I get this 1Hz frequency and 48000 points, resulting in a very "fat" logarithmic spectrum, like in the picture below.

How can I reduce resolution and thin out my plot without reducing sampling rate?

Thanks in advance.

A simple way is to split the data into smaller N sample segments (pick N to suit your desired frequency resolution, then average the power or amplitude of all the FFT spectra:

$$
  \hat X = \frac{1}{N} \sum^{N-1}_{n=0} X_n  X_n^*
$$

Where $X_n, n \in \left\{0:N-1\right\}$ is your series of shorter FFT results and $\hat X$ is the averaged power spectra.

If I remember correctly, for noise signals, the variance is reduced by $\frac{1}{N}$.


Blind source separation with one source rotated I have two images:

$$
I_1 = w_{11}A + w_{12}B, \\
I_2 = w_{21}A + w_{22} \bar{B}.
$$

$A$ and  $B$ are unknown. $\bar{B}$ is ${B}$ rotated by 180 degree. For both images, $A$ has higher signal-to-noise ratio. What would be an effective way to extract $A$ out of these two images?

Consider 2 points which are rotated version of each other, for these 2 point you could write 4 equation with 4 unknown and solve this system of linear equations (b=A*x, where b obtained from I1 and I2 at those 2 symmetric points and x are values of A and B at those points and the matrix A is obtained from Wij). Now you could solve this system of 4 linear equations for each pair of pixels, but by doing this you can't take into account the noise level unless you repeat your measurements!

To take into account the noise levels, you could use some measure of noise level like Total variation. You could build a cost function by adding the norm of (b-A*x) to a measure of noise level of your images then trying to find the unknowns by minimizing the cost function. In this case you cant use only a pair of points and you have to write the system of equations for all of points.


What is the main difference between the acoustic and electromagnetic waves? For simulation purpose, what will be the difference between using acoustic and electromagnetic waves? Just the carrier frequency only??

Reciprocity, particularly in Air acoustics.  If you are down wind from from an acoustic source like an animal and you hear it, you are less likely to be heard by it.  The same effect is true with currents in underwater acoustics but to a lesser extent because sound velocity is relatively greater than fluid velocity.

The variation of sound speed tends to be more important in simulations. There are specializations of the wave equation, actually more typically the Helmholtz equation to account for things like $c(z)$, depth dependence.

Acoustic waves do not have polarization, but do have vector particle velocity.

At the other extreme, the speed of light can not be exceeded while the speed of sound isn’t a universal constant.

The Friis equation is not used in acoustics.  Isotropic lossless propagation in acoustics is usually inverse R squared. The term “path loss” have different conventions. The Friis equation implicitly includes the aperture of the receiver, so for isotropic lossless propagation, the path loss is proportional to frequency.  Radio engineers often (stubbornly) claim that low frequencies propagate further because of the Friis equation. Lower frequencies do often propagate further but typically it is because refraction is related to wave length, so in terrestrial point to point communications, the curve of the earth, scatter by objects like houses, penetration into buildings, and mode coupling between the earth’s surface and ionosphere favor low frequency propagation.  High frequencies tend toward line of sight.  Friis isn’t why.

In acoustics, performance is typically ambient noise limited while in radio,  the limit is usually related to the antenna temperature. When an array is moving through water, flow noise can exceed ambient noise. The calculation of array gain differs accordingly.

Acoustic waves are often encountered in combination with sheer waves.  There are also slow waves in EM propagation but sheer wave propagation is essential in seismology.

In biological systems, first principle bulk properties like the Young’s modulus of living lung tissue tend to be avoided because they are very complex to model and directly measure, particularly if you don’t want to kill the subject.

There is also a distinct area of research relating acoustics to perception.  Humans (and animals) are better at localizing sounds than classic array theory predicts.


How Can Convolution and Deconvolution be Defined for 3D Images? I am trying to understand how convolution and deconvolution can be represented for 3D images/ stacks of data. I would prefer it, if you built the these concepts from 1D vectors to 3D matrices in terms of their equations. I am currently writing a report and am trying to represent 3D deconvolution as an equation.

Please do not hesitate to ask me any questions to clarify this question.

For one variable, we have

$$
y(i) = \sum_m x(i-m) \cdot h(m).
$$

For two variables it's

$$
y(i,j) = \sum_m \sum_n x(m,n) \cdot h(i-m,j-n).
$$

For three:

$$
y(i,j,k) = \sum_m \sum_n \sum_p x(m,n,p) \cdot h(i-m,j-n,k-p).
$$


Calculating Signal Approximation Error (Theoretical signal - Real hardware signal)? I've got a FIR filter, which I developed in IC. Now, when I got my results, I ended up with two signals:

Theoretical signal calculated through MATLAB as a response to my input.
Real response given by the hardware.

What is the best method of calculating the error of the real response with respect to the theroetical one?

I could go

delta = abs(Y) - abs(Yreal);
max_err = max(delta);

but I don't know what to do with the maximum error? Any suggestions?

There are lots of different error statistics you depending on what you want to know. If what you care about is the maximum error then what you suggest is fine.

Probably more common is to use the RMS error i.e. $ \sqrt{\frac{1}{N}\sum(Y-Y_{real})^2} $. This is related to the variance of the system. I'm not sure if this is exactly what you want for your case as I suspect there will be strong spatial correlation across the filter. This is something you should look at.

Another issue is that the error is dependent on the input signal. You may wish to look at the frequency response of both systems to try and remove this issue.


Why can adaptive IIR filters result in unstable solutions? For adaptive filtering, both finite and infinite impulse response (FIR/IIR) filters can be utilized. As an advantage of FIR filters in this context, guaranteed stability is often mentioned, while IIR filters do not share this property (see this related question and its answer).

I understand IIR filters are unstable when their impulse response diverges, meaning a transfer function pole lies outside the unit circle (has a radius $r>1$).

My intuition is that since in the example of system identification, only stable systems are practical to identify, it would make sense that estimates of such a systems' transfer function were stable as well.

My question is:

Why is this not the case? How can the adaption of IIR filters result in unstable solutions?

Please note that I am not asking about measures to ensure stability in IIR adaptive filtering.

The IIR filter doesn't have to be unstable, but it has the potential of being so; unlike the FIR case which doesn't have even the potential.

One reason for the (potential) unstability of an IIR (adaptive) filter is the numerical issues due to coefficient quantization. When the poles are closer to unit circle this will be critical.
This is especially important if you are using coarse quantization (like in an legacy 8 bit system?) or you are forcing the limits of your numerical precision.

Furthermore, during the adaptation process, the chaotic behaviour of the input (which is possibly nonstationary) may lead into unbounded updates of the coefficients. When the errors get quite large during a transient stage, the resulting update on the coefficients can also be quite large leading (depending on the adaptation algorithm) the filter to fall into an unstable region...