Removing transients in highpass filtering with MATLAB I've made a simple first order IIR highpass filter with a zero at z = 1 and a pole at z = 0.9. Its frequency response looks like this:

Now, I filter a DC signal using this filter. Here's the MATLAB code I use to do it:

b = [1 -1]; % Zero at z = 1
a = [1 -0.9]; %Pole at z = 0.9

figure(1)
freqz(b, a)

t = 1:100;
x(1:length(t)) = 1; % Constant function

y = filter(b, a, x);

figure(2)
plot(t, x)
xlabel('Time');
ylabel('Input Signal');

figure(3)
plot(t, y)
xlabel('Time');
ylabel('Output Signal');

As my filter is highpass, I expect the DC to become zero, or atleast become severely attenuated. However, the output signal I get looks like this:

From my understanding, this exponential output is a transient produced because I haven't set the initial conditions correctly. Sure enough, setting x[-1] = 1 solves the problem. However, this works only for this particular input DC signal. For any general input signal, how do I set the initial conditions so that transients aren't produced?

Edit : I'm aware that the filtfilt() function does forward-backward filtering with transient minimization, but I really want to port the filter to an embedded platform, so I need to understand how transient removal works. Thanks in advance for the help!

Edit 2 : As suggested by Kuba Ober, I tried setting x[-1] as the value that it actually should have been. It works fine for a DC input, but here's what happened for a sinusoidal input:

clc; clear all;
p = 0.9

a = [1 -p]
b = [1 -1]

n = 1:100; % Samples
f = 0.2; % Frequency in Hz
Fs = 10; % Sampling rate in samples per second
t = n/Fs; % Time axis

x = sin(2*pi*f*t);

% Filter with the appropriate initial conditions
y = filter(b, a, x, filtic(b, a, [], [sin(2*pi*f*0)]));

figure(1)
plot(t, x)
xlabel('Time');
ylabel('Input Signal');

figure(2)
plot(t, y)
xlabel('Time');
ylabel('Output Signal');

Here's the input signal :

And here's the output :

The first peak is visibly smaller than the second, which indicates some transients being present. I'm not entirely sure about this, but I think the reason it doesn't work is because just setting x[-1] is not enough, I also need to set y[-1]. The problem here is that there's no way to find out what y[-1] actually should be.

Edit 3 : Let me provide a little more info on the problem I'm working on. I'm trying to use filters to remove noise from ECG (Electrocardiogram) signals in an embedded platform. Here's a typical ECG signal, after filtering:

Here's what an ECG signal looks before filtering:

Note the DC offset in the signal before filtering. For filtering, I need a notch filter to remove high frequency power line noise and a highpass filter to remove the DC and the low frequency "drifting" of the signal.

The filters I use need to be linear phase, since the time domain morphology of an ECG signal is very important for diagnosis. However, my filter doesn't need to be real-time, as I'm doing the processing offline after acquiring the ECG signal from the patient. So, for implementing nonlinear phase IIR filters, I'm currently using forward-backward zero phase filtering.

One opinion that's shared by @Matt L. and @Drazick is that transients are unavoidable in real-time filtering and that I should use a longer input signal and crop off the first few seconds of the filtered output instead. This is something I'd like to avoid, as acquiring long ECGs from a living patient is somewhat difficult. Also, I do not have to filter in real-time, so any technique of transient removal that hinges on knowing the entire signal in advance is perfectly all right. Any help is appreciated!

Your first order filter recursion for some real constants $a,b,c$ is

$$ y[n] = a x[n] + b x[n-1] - c y[n-1] $$

with the two initial memory states $x[-1]$ and $y[-1]$ at $n=0$.

Your "no transient" condition can be translated to $y[0]=0$ and a necessary second condition so that you can solve for both of your memory states. The second condition could be, that the discrete derivative of $y$ also vanishes at $n=0$, so $y[0]-y[-1]=0$. You can also take any other condition that seems sensible to you.

The two equations give you a unique solution for the two unknown memories, namely:
$$y[-1] = 0$$
 and
$$x[-1] = -\frac{a}{b}x[0] $$

Alternatively, your conditions may better be chosen as $y[0]=0$ and $y[0]-y[-1]=x[0]-x[-1]$ in order to capture the initial slope of the input. The resulting recursion equation at $n=0$ is then
$$0=a x[0]+b x[-1]+c(x[0]-x[-1])$$
giving you the solution
$$x[-1]=-\frac{a+c}{b-c}x[0]$$
and
$$ y[-1]= -\left(1+\frac{a+c}{b-c}\right)x[0]$$

(Please check my calculations!)

But in general you cannot expect a simple initial condition to give you the same result as knowing the signal history. So you can only take this to a certain point and in general it would probably be better if you just discarded the transient response of your output.


confusion sampling vs quantization? while converting analog signal to digital equivalent,we have a process that is called analog to digital conversion and it has two main steps/stages sampling and quantization? I am confused whether y axis is discretized in sampling stage along with x axis or only y axis is dicretized in quantization stage since x axis is already discretized in sampling stage?

Initially i thought that sampling only involves discretization of x axis but now i was watching following youtube lecture where it tell at approximately between 10 min and 11 min that discretization of both x and y axes are done during sampling stage?

https://www.youtube.com/watch?v=xUCsfKA8bi0&list=PLm_MSClsnwm9I2iviE0YKt6PZTyQCYc8j

Sampling is the process of making the x-axis (time) discrete and quantization is the process of making the y-axis (magnitude) discrete.  You can sample without quantization (such as done with an analog sample and hold circuit). Quantization is introduced through rounding or truncation when the sampled analog signal is mapped to a digital representation.

Sampling only causes aliasing distortion to the signal; if there is no energy above the half the sampling rate for real signals with real sampling, or above the sampling rate for complex signals with complex sampling, absolutely no distortion occurs assuming a perfect sampling clock with no hitter/phase noise (both aspects are not physically possible to totally eliminate but can each be reduced to below our level of concern).

Many A/D's, such as successive approximation converters, are implemented with a fast analog sample and hold so that the sampled signal is held while the conversion to digital takes place.

Quantization noise adds noise to the resulting signal, well approximated in most conditions as a uniform white noise: the noise distribution in magnitude is uniformly distributed over the magnitude of a quantization level and is white in frequency to the extent that each sample of the noise is independent of the next.

From that approximation of the noise we can arrive at the SNR estimate for a full scale sine wave to be 6.02 dB/bit + 1.76 dB, which we can increase further by “over sampling”; since the noise added is white, we can filter out the noise in the band we don’t need.

I have detailed that process further at this post: What are advantages of having higher sampling rate of a signal?


Definition of average power? There are two kind of average power I encountered in random signal class and textbook:

definition 1: average power =$$E[|x(t)|^2]=R_{xx}(0)=\int^\infty_{-\infty} S_{xx}(f)\,df$$
definition 2: average power =$$\lim_{T\to\infty}\frac{1}{2T}\int^{T}_{-T}|x(t)|^2\,dt $$
It seems that the two definition is different unless the signal is ergodic in $x(t)^2$. By the way, our teacher keep using the term "assuming average power of white noise $n(t)$ is $\sigma_N^2$" in test. But from the above definition, the average power is infinite for continuous white noise. I think in discrete white noise the average power can be well-defined but in continuous one I am wondering our teacher may take the definition wrong.

The first definition works for deterministic as well as for random signals. For random signals we define the autocorrelation by

$$R_x(\tau)=E\{x^*(t)x(t+\tau)\}\tag{1}$$

where $E\{\cdot\}$ is the expectation operator. For deterministic power signals (i.e., signals with finite but non-zero power and infinite energy) we define

$$R_x(\tau)=\lim_{T\to\infty}\frac{1}{2T}\int^{T}_{-T}x^*(t)x(t+\tau)dt\tag{2}$$

From $(2)$ it is clear that for deterministic power signals both definitions of average power are equivalent, because for random signals as well as for deterministic power signals we have

$$S_x(f)=\mathcal{F}\{R_x(\tau)\}\tag{3}$$

Also take a look at this related answer.

You are right about the power of white noise. In continuous time, white noise has infinite power. Only discrete-time white noise has finite power. In sum, teachers and professors are only humans.


what is better: up- or downsampling? I want to compare two signals or curves.

Unfortunately they have different sample rates. The first one is sampled at 30 Hz, the second at 2000 Hz. Matlab has the function 'resample' and I  thought, this would make comparison much easier.

My question is: Is it wiser to downsample the second curve or to upsample the first one?

Edit:
I did as I was told. In the picture on the left are the original curves. The curves on the right are resampled. Top right is upsampled, bottom right is downsampled.

I know the curves are of different length, when resampled. However, the second signal will get cropped to the first signals length.

In short:

Upsampling: does/should not loose information  (if done wisely), then safer,
Downsampling: may loose information (if done unwisely), yet more computationally  efficient.

So if you compare data at different rates, and in an evaluation phase when one tries to define how the comparison should be done (which features are compared, with what metric, with which externals, like precision, efficiently, robustness, etc), a very basic first approach would be to upsample both signals by an integer factor to reconcile their sampling. Here, this is quite OK, 6000 Hz could work for both.

The reasons could be, for DSP practice:

with integer upsampling, you are not obliged to use tricky filtering techniques, and simple linear interpolation is easy. Remember though that upsampling adds some information.
with signals are the same scale, you can extract features and compare them: noise, variability, slopes, derivatives, etc. You can apply the correct scale/offset correction (they don't have the same amplitude), use linear or non-linear transformations (Fourier, etc.) to check whether similarities features exist in different domains
at this stage, you should be able to quality and quantify the features/metrics needed to build your algorithm. And you can estimate what you can loose from both data. For instance, you can create your algorithm, and see how it is robust to signal downsampling, and how far you can go. Here, your algorithm begin to be more efficient.

Once you have spent some time on that, you are more ready to start from scratch again, and decide whether you should downsample or not, with which method, etc. to reach both safety and efficienty.

To provide alternatives, if you decide to compare data on the Fourier domain, you can easily cast two signals with different length and sampling rate to the same number of Fourier coefficients, helping a comparison in the frequency domain.


Understanding PSD: Why Does Power at High Frequencies Affect Low Frequencies? I'm trying to wrap my head around power spectral density on a conceptual level, but I am having some difficulty. Suppose I have a communication system where I am receiving and sampling white Gaussian noise, which has a uniform PSD. Increasing the bandwidth of my receiver results in greater noise power across all frequency components. I can follow the math of how this occurs: integrate the PSD across more frequencies and you will increase the total power, which is uniform for white noise. However, how does it conceptually make sense that letting in higher frequency power components directly increases the power of the low frequency components?

EDIT: Some commentators have pointed out that my question was worded poorly. I apologize. I will try to clarify what I mean by turning it into a concrete thought experiment. Suppose I have a LPF with bandwidth $B_1$. Then, the variance of the DC component in the DFT domain after sampling WGN will be some $\sigma_1^2$. Now, compare that with a wider LPF with bandwidth $B_2$, where $B_2 > B_1$. Then the DC component variance would be greater, $\sigma_2^2  > \sigma_1^2$. Now, what if instead I implemented this second LPF filter a different way, with a parallelized structure: the sum of 2 different filter paths. The first filter path is the LPF with bandwidth $B_1$ and the second is a bandpass filter with passband $[B_1, B_2]$. This is followed by sampling. This is an equivalent structure to the LPF with bandwidth $B_2$, so the output statistics should be the same. However, notice that the second path is DC blocking, so it will never have a DC component. This should mean that the variance of the DC of the sum of the paths is equal to just the variance of LPF path, $\sigma_1^2$. But clearly, despite this being an equivalent filter structure, the output statistics are not the same. Why?

The spectrum allows you giving another complementary look on the data.
It is the other side of the same coin.

In Frequency Domain we are looking on the energy of each spectral component of our data.
Each component is independent of other component in the meaning each component is orthogonal to any other component.

In the case of White Noise, the PSD is the Fourier transform of the Auto Correlation function.
The model is Continuous White Noise is one which has infinite power.
Its power becomes finite once it is sampled by a sampling chain which has Low Pass Filter (Which the component which limits the incoming power) and AD.

Now, just like LPF will limit the power of any signal which has bandwidth which is higher the LPF bandwidth it does so to White Noise as well.

Namely, the Constant Level of the PSD is a property of the White Noise alone (In case of Gaussian White Noise it is the $ \sigma $ parameter) while the power of the sampled noise is derived by the level of the PSD and its bandwidth (Which is a property of the LPF).

There is no high frequencies which changes the values of low frequencies.
The power gets bigger because of the hidden assumption that when you sample with higher rate you use LPF with higher bandwidth which lets more frequencies in and hence more power.


Memory efficient filtering with scipy.signal in Python I have several EEG channels that need filtering, and one channel may reach a recording length of 12 hours or potentially even more. This results in billions of data points stored which, when loaded, is represented as a NumPy array (by default as 16 bit long integer). In order not to load all data into memory at once, I make use of Python's mmap. This solution works until I call scipy.sosfiltfilt(). As I observed, this method creates copies of the entire signal which gets eventually filtered and is returned. The problem arises with really large arrays when making copies becomes simply unfeasible in terms of memory usage.
I have done some research, however, I have not come across a solution. Therefore, my question is whether there is any solution to filter such lengthy signals efficiently. My goal is to achieve this by for example, using an in place updating approach instead of creating copies without having to reimplement the method. Or does scipy.signal support mmap in any way? Or should I just create smaller chunks of the data? Or is there any other feasible way to do this?

The way to do this is to break the signal up into chunks and process each one at a time.
I asked in a comment if you could accept one-dimensional filtering, but I guess you can do bidirectional filtering, too; you just have to do it yourself.
For one-directional filtering:

Break up the signal into chunks that can fit in memory  (Is there a convenience function for this?  Maybe np.nditer?)
zi is initially all zeros
For each chunk:

Run it through scipy.signal.sosfilt, with zi=zi
Set zi = zf (from the output of the filter)
Write the chunk of output filtered data to your new array

When you've assembled all the output chunks, it will be as if you had done one continuous sosfilt in one direction.
For bi-directional filtering:

Do the whole thing again, but in the reverse direction.

For this to be exactly equal to sosfiltfilt, you probably need to keep some additional tail past the end of your signal, and then use that as your initial conditions when going backwards?  Look at the code for sosfiltfilt or filtfilt and confirm how they work.
If you are ok with FIR filtering, then that's easier: You can make a zero-phase filter and just do that in one direction (using lfilter instead of sosfilt), with the same process as above, and then throw away the head of the output, so that it looks as if it were done with no offset.


Finding PHASE between two sinusoids with FFT/DFT (ArduinoFFT w/ C) I am having trouble finding a known phase between two data sets, lets call them data set 11 and data set 12.
I am obtaining a sample size of 1024 samples, at a sampling rate of 500kHz with a total sample time of 2.048ms.
All data and code used in this question can be found at this repository - https://github.com/smake5730/fft_phase/graphs/contributors
As can be seen from the sample data below the two data sets show the voltage is leading the current by around 90°. (zoomed in)

When I compute the FFT with 1024 bins I get clear peaks. (zoomed in)

I apply the Hanning window during this and do not use zero-padding.
I use the location of the largest peak with the atan2(a/b) method to obtain the phase and find that I get the following below. The first data gives me +90° and the second data gives me -90°. The frequencies are correctly identified.

I have tried adding zero-padding to the end of the signal and it does not solve the issue. I have also tried upsampling and it does not give the correct phase either.
Does anyone have any ideas/hints on what I could try or where I am going wrong?

EDIT
As was pointed out below, my code takes individual bins and incorrectly matched bins were used in the second set.

You can't do better than this for accuracy in my experience.

A Two Bin Solution

Select a frame size of 1 1/2 cycles or 2 1/2 cycles.
You only need to calculate two bins, either 1 and 2, or 2 and 3.  No need for a full DFT.  So it is very efficient as well.  By placing the frequency squarely near the center of the two bins you minimize noise effects and avoid needing to deal with near bin frequencies which call for a different set of equations for the fitting vector.
The alternative form is for when the frequency is really close to a bin.

DFT Bin Value Formulas for Pure Real Tones
An Alternative Form of the Pure Real Tone DFT Bin Value Formula

There are other tricks to do it point by point.  I'll let others provide those.


How to represent multiple related signals via mathematical function? I have three data sets from same source. I know they are interdependent and stochastic also. The datasets are of very high frequency (see the graph, X axis is time and Y axis is value). I want to represent these three signals in a single mathematical function. Is it possible?

In the above graph, please see there is a pattern or seasonality. If we change one then other also changes. I have tried to google for it but did not find any model that can completely represent this type of data in a single expression.

Try the Eureqa software. It tries to find the simplest mathematical relation between multiple data. If the relation truly has a simple symbolic form, I think it should find it. Otherwise it will produce some strangely formulated approximation.


Most efficient way to extract frequency components of audio signal with near-human accuracy in real time I'm trying to work out how (if possible) to extract the frequency components of an arbitrary audio sample (typically music) in an FFT-like manner, but in my research on the FFT algorithm, I'm learning that it suffers some severe restrictions for this purpose.

There are 3 problems which the FFT is presenting:

Because the FFT bin-resolution is equivalent to your window size, to achieve a pretty reasonable accuracy (say 1 Hz), you need an unreasonably long window (say 1 second). This means you can't detect transients or newly introduced frequencies quickly. It also means that the problem can't be solved with a faster CPU and a higher sample-rate -- the restriction is intrinsically tied to time.
Humans perceive frequency logarithmically, but FFT bins are spaced linearly. Eg a difference of 20hz at the low end of our hearing is huge, whereas a difference of 20hz at the high end is imperceptible. So to get the accuracy we require at low frequencies, we have to calculate far more than we require at high frequencies.
Some of these problems can be resolved by interpolating between FFT bins. This may work for much musical audio, because frequencies will often be spaced quite far apart and so no more than 1 frequency will leak into a pair of bins. But this will not always be the case, particularly for inharmonic sounds like percussive instruments. So interpolation is really just guesswork.

From what I understand of the DFT/FFT algorithm, the outputs (bin amplitudes) are effectively the correlation of the sine/cosine at each bin's frequency. It strikes me that if the algorithm could be redesigned so that the bin frequencies are spaced non-linearly (i.e. we correlate a different set of sines/cosines), then we could achieve pyschoacoustically-equal resolution at all frequencies. Is this possible, or is it a pipe-dream based on my incomplete understanding of the maths involved?

I guess I could also solve the problem with brute-force, by correlating sines/cosines at every single frequency I'm interested in. I'm not too clued up on the maths here. Is this possible? What sort of efficiency? Would it solve my problem?

Is there a different way to achieve more accurate, real-time, frequency decomposition of a signal? CPU efficiency is a concern, but not the major concern -- I'm partly interested in whether it can theoretically be done at all. However, something that's feasible in realtime on a modern desktop machine would be ideal.

First, with the classic short-term Fourier transform approach, there are alternative to interpolation - in particular techniques making use of phase information to recover the instantaneous frequency (See this question) which can give you very accurately the position of a spectral peak without an increase of FFT size. The drawback, as you correctly said, is that you are not increasing the ability of the system to discriminate adjacent peaks - but this is already a great improvement compared to using the central frequency of the FFT bin index.

Regarding your brute force approach... FFT has this resolution limitation (with bins at $\frac{sr}{FFT\_size}$) because it is "probing" the signal with complex exponentials that have a null integral over the analysis window - and these must have a period which is an integer division of the length of the analysis window. If you try to naively use the same approach to "probe" frequencies at other ratios, it won't work because you'll be "probing" your signal with functions that do not contain a complete cycle. If you wanted to make it work, you would have to increase your analysis window to the least common multiple of the periods of all the signals in your brute-force search, and this is precisely what you want to avoid!

There is another brute-force approach that works: "probe" your signals with windowed complex exponential (Gabor wavelets). These are characterized by a center frequency, a center time, and a bandwidth (which measures how the wavelet is spread over time or over frequency). You'll have to evaluate many, many, many correlations between your signal and these wavelets at as many time offsets, frequencies, and bandwidths you want. The result will be the same as a very flexible "tiled" STFT in which an optimal window size is selected for each time-range and each frequency-band. Besides the computational cost, the downside is that there is no efficient algorithm, and no causal algorithm (you will need to know as many samples in advance as the longest wavelet in your dictionary). If you want to experiment with these techniques, have a look at MPTK.

A last class of approaches working quite well for music signals are subspace tracking / high resolution methods. They work the problem from a model estimation angle - "assuming this signal is a sum of $k$ exponentially damped sinusoids, what are the best amplitudes/frequencies giving the least square error?", and as such the number of samples required to get an accurate estimate is only dependent on the noise level. There are some hurdles to make them work, though:

The number of components $k$ must ideally be known in advance. Order estimation methods like ESTER have been proposed to optimally pick the model order.
They perform well in the presence of white noise - this requires the signal to be whitened before analysis; performing the analysis in individual channels of a filter bank helps too.

These are computationally expensive, but they can work online, with short windows if model orders and/or noise are low.


How to Zero Pad in Order to Perform Filtering in the Fourier (Frequency) Domain? Consider an $M\times N$ image $f$ and an $G \times K$ filter $h$. Given that convolution in the spatial domain corresponds to multiplication in the Fourier domain, then we can perform a convolution of $f$ with $h$ (that is, we can filter image $f$ using filter $h$) by multiplying the Fourier transforms of $f$ and $h$, $\hat{f}* \hat{h}$, and then perform an inverse discrete Fourier transform of the result, $\hat{f}* \hat{h}$, where $\hat{f}$ and $\hat{h}$ are respectively the Fourier transform of image $f$ and kernel $h$.

However, to do that, we need to make sure that $\hat{f}$ and $\hat{h}$ have the same dimension. So, we use the zero padding technique. We pad both the image $f$ and $h$, before finding their Fourier transforms. More precisely, we pad each of them with zeros, such that their new size is $(M + G - 1) \times (N + K - 1)$. This padding is apparently required (for both $f$ and $h$), but I am trying to understand why.

Why exactly is this the case? I think it has something to do with the assumption that $f$ is periodic. So, the conclusion would be that we can't simply zero pad $h$ (without also zero padding $f$), but I don't get why.

The result of a convolution of a data vector of length M with a kernel of length G is of length M + G - 1.  (the maximum length of the non-zero portion, even though the limits of integration is sometimes written as from -infinity to +infinity)

This is clearly (G - 1) elements longer than the original data vector.

So where do these new, "extra", additional result values go?  If you don't zero pad when doing an FFT fast convolution, these "extra" result values simply get added into some portion of the first part of the length M convolution result, wrapping around the edge, which likely messes up your desired linear convolution result.

But if you zero-pad, then these "extra" values get added to the zero-padding, which is less likely to wrap-around far enough to mess up your desired convolution result.

Both are circular convolution, but a circular wrap into a bunch of zero-padding is far far easier to delete or remove, then when it's added into, and thus all mixed up with the results you most likely actually want.  e.g. becoming artifacts along the edge of your filtered image.

You have to zero pad the kernel to the same length as the data so that their FFTs also end up the same length, which allows straight forward multiplication of those 2 equal length FFT vectors, before IFFT-ing the multiplied result.

So, how much to zero-pad?  By any amount >= (G - 1).  Perhaps rounding up to the nearest size for which your chosen FFT implementation is performant, usually a size which is the product of only very small primes, such as 2 (and/or maybe 3, 5 or perhaps 7 as well).


Peak detection approach What are the peak detection algorithms in existence? I have noisy data, and I would like to implement peak detection for this data. Data is in reverse, actually I am trying to determine the bottom.

Here is a snapshot of the data in Excel. I like to detect both bottoms. I thought about passing the data through a low-pass filter and subsequently do a moving average where I determine the peaks and within the moving average I do another search. I have zero DSP background; this is just a common sense approach. I would like to hear what the experts recommend.

Ktuncer, there are a number of methods you can use here. One method that I would recommend is to use a Discrete Wavelet Transform, (DWT), and in particular, look at the Daubechies Wavelet. I would pick, say, Daub-14 / Daub-Tetra.

Basically what you really need to do is 'trend' your signal, and then from there, do a min or max pick. That will get rid of your outliers. A daub-14/daub-tetra wavelet transform can help you do this, and this helps especially since you do not know the nature of your signal. (Using daub-14, you can accurately represent polynomial signals of degree 14/2 = 7, and it looks like you wont need more than that).

The computation of this wavelet transform essentially 'compresses' your energy into a few indicies. Those indicies represent weights on basis vectors. The rest of the weights will (ideally) be near zero. When you have noise in your signal (as you do), those weights that were normally zero have some weights now, but you can simply zero them out and 'denoise' your signal. Once that is complete, you can then do a simple max/min detect.

There are more details involved, you may send me an email if you like to discuss how to implement it. I have done similar work on this before.

EDIT: Here are some images illustrating Daub-Tetra Denoiser:


Looking for low-noise audio data to benchmark some noise-processing techniques I'm looking for a source of very low-noise audio data for use in benchmarking various methods of random noise correction methods.

My intention is to have a very clean signal to which I can add known random noise such that I can test the efficacy of various methods of dealing with that noise.

It would need to be in some lossless format (preferably .wav).

Does anyone have suggestions?

Thanks!

That is a simple question:

All sound that passes into analogue domain using digital audio converters (DACs) have a SNR inferior to about 113 at best, and more probably 96-100... everything recorded using microphones will have inherant noise, unless the recording has been copy pasted onto a silent track leaving intervals of silence.

So what you need to do to have a very clean signal, is to take a totally digital sound which was generated with i.e. sine waves and has not had reverb and echo effects that raise the noise in the silent time(s) of the recording, and digital synthesizers which return to 0 after the note envelope. You can go though a track that contains elements of silence in Audacity or a more easy to zoom editor, and check that you are having times in the track where the amplitude = 0.0; or program a fast test code to count that the track has amplitudes of 0.0 lasting for more than 500ms or something.

if you take any recording and vary the volume of the track, for example you take an AC-DC rock song, and apply an artificial envelope that fades the track to 0.00 regularly i.e. every couple of seconds, it also makes a test file with totally silent intervals.

you can also write a code that sais : if volume < than 0.05, volume =0; to artifically make any track you want totally silent at times.

find digital recordings on freesound and measure them for 0.00 intervals in between sounds, check digital music artists that boast about 192khz recording equipment for a particular album, get a feel of which artists use state of the art production equipment, check standard commercial tracks that are used to measure monitor speakers, there are many options.


Channel coder: Does encoder help? Input signal distortion: I am studying the channel coding theory. The channel coder adds redundant information to the input signal, that helps a decoder to reconstruct the signal.
BUT: Can encoder help to eliminate effects such as Doppler, offset...?
Edit 1
when we are talking about improving encoder in practice, what does it mean? What should be improved?
Edit 2:
Why has an interleaver complicated implementation in practice? How to reduce complexity?
Edit 3
S/P or P/S converters and interleaver, do they make the same task?

BUT: Can encoder help to eliminate effects such as Doppler, offset...?

No, since the encoder doesn't "see" them. The overall system, including en- and decoder can be optimized to deal with such effects.
But honestly, something that corrects a Doppler offset would be called "frequency synchronization", and work on the digitized signal, not on codewords, so it's not much to do with channel coding.
There is techniques that you find both in decoders and for example equalizers, but those are really more shared methods, not "channel coding". Rule of thumb would be for me: if it's something in a receiver that has to do with optimizing the signal prior to mapping it to information, it's digital signal processing, afterwards, it might be channel or source coding.

when we are talking about improving encoder in practice, what does it mean? What should be improved?

This is a  bit to broad, because nothing is ever optimal for all use cases. Example: One man's block code is "optimal" if it's as large as possible to come as close as possible to theoretical limits in code gain. Another man needs to decode something within a very short time, and can't wait for a gigantic code word to be completely received. These two people will have very different ideas what an "optimal" thing is.

Why has an interleaver complicated implementation in practice? How to reduce complexity?

Interleavers are typically among the least complex things in a receiver or transmitter, so not really sure what you're referring to.
NB: there's often things that can be tackled with knowledge that the decoder has; any decision-guided synchronization algorithm is testament to that. In some cases, simple types of ISI might be corrected by a convolutional decoder - but that's really, like, a special case that you can't generally rely on.


DFT (FFT) of a Real Even Function Doesn't Yield Real Only DFT Signal There are a lot of similar questions already, but I can't find one answering this problem.

I used the following MATLAB code to perform the FFT on a Lorentzian input signal

x = linspace(-10,10,2^10);
y = 1./(x.^2 + 1);
F = fft(y);
IF = imag(F).^2;
power = trapz(IF); %integrating the imaginary part squared

Now, from everything I've read, the FFT of a real, even function is supposed to give $0$ imaginary part. But in this case, I get

power = 3.86

which doesn't seem to be just a small numerical inaccuracy. What's going on here?

Attached is the imaginary spectrum after fftshift.

When you work with DFT you need to remember one big assumption of the DFT - The samples are part of periodic signal.

For example, assume your signal in on the grid vX = [-3:3].
Let's say it is the simplest symmetric function: vY = vX .^ 2.
If you apply fft on it the result won't be real.

Why is that?
Because the DFT assumes the input is periodic.
Hence the input signal is something like [vY, vY, vY, ..., vY].

Moreover, for the DFT the indexing is [0, 1, ..., N - 1] which in this case is [0:6]. Now if we do periodic continuation, what will you get on the grid of [-7:6]?

figure();
plot([-7:6], [vY, vY]);

It won't be a be symmetric around 0 (Of course look at [-6, 6]).

Let's define vYY = vY(1:6).
If you apply fft(vYY) the result is pure real.
Why is that?
Because for the DFT it sits on the grid [0:5].
If you plot a periodic continuation of it:

figure();
plot([-6:5], [vYY, vYY]);

Look on [-5, 5], now it is symmetric.
Hence its DFT is pure real.

MATLAB Code

% The Grid
vX = [-3:3];
vY = vX .^ 2; %<! Even Function on the Grid

figure();
plot(vX, vY);
title('Figure 001');

vYDft = fft(vY);
norm(imag(vYDft), 'inf')

norm(imag(fft(vY(1:6))), 'inf')

figure();
plot([-7:6], [vY, vY])
title('Figure 002');

figure();
plot([-6:5], [vY(1:6), vY(1:6)]);
title('Figure 003');

Output:

  ans = 7.2678
  ans = 0

Figure 001

Figure 002

Figure 003


GMSK CFO estimation and timing recovery My goal is to estimate the CFO and timing offset in GMSK modulated signal. I synthesize samples in Matlab and my goal is to add CFO and time offset to learn how to deal with mentioned distortions and successfully recover transmitted signal. I would like to use non-data-aided methods and rely only on the signal structure. To be honest, I read a lot and I have few methods in my mind, like: M&M (but I guess its for timing recovery), squaring method or FFT method (to calculate CFO), PLL/Costas Loop (for fine CFO estimation). I don't know if they really fits for GMSK so my question is, how to tackle this problem? What should be done first in the recovery algorithm to stay in the right order? My goal is to deep into signal recovery and to learn different approaches.
Thanks in advance

For optimum coherent detection in an AWGN channel, GMSK with partial response signaling (intentional overlap of pulses to increase the amount of data within a given time and therefore bandwidth efficiency) can be done with maximum likelihood sequence detection, combined with Viterbi decoding which basically compares the correlation of all possible phase trajectories (the number of which is increased substantially due to the pulse overlap) and selects with soft decision detection the most likely set of data bits that would have created that trajectory. With full response signaling there are just two possible trajectories since only one data bit is involved, but with partial response signaling there are multiple data bits to create each trajectory due to the longer duration from each data bit that is longer than a symbol period. For example if each transition lasted 2 symbols, there would be four possible phase trajectories corresponding to the 00, 01, 10, and 11 data patterns. The trajectories are unique to the patterns and thus by choosing the most likely trajectory (after removing phase, frequency and time offsets) the data can be demodulated. As the overlap increases (duration of the phase transition over multiple symbols) the complexity of the optimum receiver significantly increases due to the number of phase trajectories and it can be seen that a trellis diagram results where the Viterbi algorithm can be used to efficiently determine the most likely path through the trellis. Further Laurent Decomposition can be used to decompose the signal into PAM modulations and reduce the possible combinations, simplifying the receiver. Other sub-optimum receiver implementations exist as well.
Here are related links for GMSK demodulation, including further details for carrier recovery, symbol synchronization and using Viterbi decoding and the Laurent Decomposition:
https://www.mathworks.com/help/comm/ref/comm.gmskdemodulator-system-object.html
https://core.ac.uk/download/pdf/195023475.pdf
https://www.atlantis-press.com/article/125913177.pdf
https://etd.lib.metu.edu.tr/upload/1057348/index.pdf


Am I understanding the phase angle from FFT correctly? Long story short, I have a sinusoidal wave with DC offset going in with a given frequency of say 100 Hz and $\phi = 0$ (which I'm guessing is really $\pi/2$ since it is a sine wave not cosine). I get an output of the same frequency and, using FFT, I would like to find the phase between these two. What I did is I found the peak frequencies (-100 and 100 Hz), found their bin and from the same bin I extracted the amplitude and the phase.

Is that correct?

If you're interested in the time delay between the two signals, it's probably most useful to compute the cross-correlation function of the two signals and find its maximum. The location of the maximum will give you the time lag.


What is the difference between Image Inpainting and Image Denoising? Both of them seem to do the same task. Wherever the pixel values within the image refers to something that was not there in the original image, we call that pixel to be affected by noise. If this noise is less, we apply noise reduction. When the pixel refers to totally random value, we apply image inpainting.

So is image inpainting just an extreme case of image denoising?

Also if a large portion of pixels get affected, we call that a "noise" and need to apply denoising algorithms; whereas when a small portion gets affected like a scratch on the image, we apply inpainting techniques. But my question is: that scratch can also be regarded as a noise only. So why dont we call that process image denoising?

Inpainting solves a binary problem of missing data. These masked or missing pixels to not contain information that is used in the reconstruction. Also, is it assumed that non missing pixels are to remain as is.

Denoising usually refers to the situation where all the pixels contain some noise component, which often also has some global prior independent of the image content. The algorithm tries to remove the noise from all the pixels, and uses the information inside these noisy pixels for the reconstruction and estimation.

The problems are not unrelated and could be reduced to private extreme cases of each other but in general they solve different types of problems.


How to determine if the system is linear or nonlinear Given a signal

$$y(n)=x(n)\cos(wn)$$

How to determine if the given signal is linear or non linear?

I am getting so confused reading the theory. If someone could help step by step to determine this. It would be great thankyou

As Matt L. says you'll need to check for homogeneity and, possibly, additivity.

Homogeneity

That test says that if:
$$
y[n] = f(x[n])
$$
then
$$
A \cdot y[n] = f(A \cdot x[n])
$$
for all scalar $A$.

Additivity

This test says that if
$$
y_1[n] = f(x_1[n])
$$
and
$$
y_2[n] = f(x_2[n])
$$
then
$$
y_{\rm tot}[n] = f(x_1[n] + x_2[n]) = y_1[n] + y_2[n]
$$

You just need to apply these for your system.

I say "possibly" because if the system fails one of these tests, then it will not be linear --- so there is little point testing the second condition.


Discrete Kalman filter for a continuous system The question is related to the implementation of a discrete kalman filter given a description of the system model in continuous time.

I will give an example. Suppose we have a mass, spring and damper system as below:

The differential equation describing this system is:

$$\ddot{x} = -\frac{k}{m}x - \frac{b}{m}\dot{x} + \frac{1}{m}F$$

Therefore, if the system states are:

$$\begin{bmatrix}
x \\
\dot{x}
\end{bmatrix}$$

Then, the state space matrices are (assuming we can only measure position):

$$A = \begin{bmatrix}
0 & 1 \\
-\frac{k}{m} & -\frac{b}{m}
\end{bmatrix}$$

$$B = \begin{bmatrix}
0 \\
\frac{1}{m}
\end{bmatrix}$$

$$C = \begin{bmatrix}
1 & 0\\
\end{bmatrix}$$

$$D = \begin{bmatrix}
0
\end{bmatrix}$$

I would like to implement a discrete kalman filter for this system, following the regular discrete system convention:

$$\hat{x}_{k} = A\hat{x}_{k-1} + Bu_{k} + w_{k}$$

$$z_{k} = C\hat{x}_{k} + v_{k}$$

Of course, the A, B, C and D matrices of the continuous system are not the same ones as for the discrete system.

I have tried to convert the continuous system to a discrete one (using ZOH, see also here) but it didn't seem to work.

Any ideas on how to tackle such a problem? An example would also be nice.

The mass-spring-damper combination is an LTI system described by the following continuous-time linear differential equation

$$ \ddot{x} = - \frac{k}{m} x - \frac{b}{m} \dot{x} + \frac{1}{m} u $$

where $u$ is the deterministic input force (N),  $k$ is the spring constant (N/m) and $b$ is the damping coefficient (N.s/m).

Assuming two states as the position $x(t)$ and the velocity $\dot{x}(t)$; $X = [x , \dot{x}]^T$, and assuming a single measurement as the position $x(t)$ and assuming no process noise, and scalar deterministic input $u$ only added to the process equation, then the continuous-time linear dynamic system process and measurement equations for the Kalman filter can be written as

$$ \dot{X} = F X + G u $$
$$ y = H X + v $$

where
$$
F = \begin{bmatrix}
0 & 1 \\
-\frac{k}{m} & -\frac{b}{m}
\end{bmatrix}
$$

is the system dynamics matrix,

$$
G = \begin{bmatrix}
0 \\
1/m
\end{bmatrix}
$$
is the deterministic input matrix,

$$
H = \begin{bmatrix}
1  & 0
\end{bmatrix}
$$
is the measurement matrix, $v(t)$ is the measurement noise.

The first step in the derivation of the discrete Kalman filter is the discretization of the associated matrices by first finding the fundamental matrix $\Phi(t)$ which is given by :

$$ \Phi(t) = \mathcal{L}^{-1} \{ (s I - F)^{-1} \} $$

where $s$ is the complex variable of the Laplace transform.

For the given $F$, the fundamental matrix is found to be:

$$
\Phi(t) = e^{\frac{-b}{2m}t}
\begin{bmatrix}
\cos(\omega t) + \frac{b}{2m\omega } \sin(\omega t)  & \sin(\omega t) / \omega  \\
-\frac{k}{m\omega } \sin(\omega t) & \cos(\omega t) - \frac{b}{2m\omega } \sin(\omega t)
\end{bmatrix}\\
$$

where $\omega = \frac{\sqrt{4mk-b^2}}{2m}$.

Then the discrete fundamental matrix of the Kalman filter is found from $\Phi_k = \Phi(T)$ where $T$ is the sampling period.

$$
\Phi_k = e^{\frac{-b}{2m}T}
\begin{bmatrix}
\cos(\omega T) + \frac{b}{2m\omega } \sin(\omega T)  & \sin(\omega T) / \omega  \\
-\frac{k}{m\omega } \sin(\omega T) & \cos(\omega T) - \frac{b}{2m\omega } \sin(\omega T)
\end{bmatrix}\\
$$

We should also find the discrete versions of $H$ and $G$ which are found as $H_k = H = [1 , 0]$ and

$$G_k = \int_0^T \Phi(t) G dt .$$

By the help of a symbolic integration package (such as provided by MATLAB) one can compute the discrete deterministic input matrix $G_k$ as

$$
G_k = \frac{e^{aT}}{m (a^2+w^2)}
\begin{bmatrix}
\frac{a}{\omega} \sin(\omega T)-\cos(\omega T) + e^{-aT}  \\
 a (1-\omega) (\cos(\omega T)-1) + (a^2+ \omega) \sin(\omega T)
\end{bmatrix}\\
$$

where  $ a = -b/2m  $ and $\omega = \sqrt{4km-b^2}/2m$ as before.

After having found all the necessary discrete matrices $\Phi_k$, $G_k$, $H_k$ we can setup the discrete Kalman process and measurement equations through which the recursive solution can be applied.

$$X[n] = \Phi_k X[n-1] + G_k u[n]$$
$$ y[n] = H_k X[n] + v[n]$$

where $u[n]$,$v[n]$ and $y[n]$ are sampled versions of deterministic input, measurement noise and measurement respectively, and $X[n]$ is the sampled version of the state vector.

Then you will apply the Kalman recursion to estimate the states. Below are some simulation results from a MATLAB implementation of the associated Kalman filter. In the simulation, the true signal is generated by solving the associated ODE with a simple Runge-Kutta numerical method. Eventhough it's not required for this LTI system whose analytical solution is already available, it'll be helpful to use a numerical method later if more complicated (nonlinear) differential equations are encountered.

In this simulation, the model was exact. And therefore it's a little trivial. For more realistic applications, one should either replace the linear differential equation with more realistic nonlinear ones hence change the model and/or apply process noise to compansate for the lack of exact process model information.


Transforming digitized noisy signal before applying cross-correlation I'm trying to grasp the concept of cross-correlation as it applies to CDMA in the GPS C/A signal where noise is involved and the SNR is low.

My understanding is that before calculating the cross-correlation with the C/A code one has to transform (just shift by its average value?) the digitized and demodulated received signal so that there are both positive and negative values in the sequence and the sum of these on average should be zero, is this correct?

(I understand that there's much more to acquiring the GPS C/A signal, but I'm trying to grasp how cross-correlation can be applied to a digitized signal)

To offset or not is simply in your definitions on how you want to do the math involved and can be convenient for further processing.  Subtracting a constant, or otherwise scaling the waveform, does not change the signal to noise ratio.

Correlation is to multiply and accumulate, and the cross-correlation and auto-correlation functions show this correlation as the sequences are shifted in time relative to each other.  Cross-correlation is between two different functions, and auto-correlation is between a function and itself.

The cross-correlation of x[n] and y[n], with a shift index m, can be given as:

$$\rho_{xy}[m] = \Sigma_n x[n]y[n+m]$$

(Note: For complex signals, the product must be complex conjugate product)

With that said, consider this very simple example using an 11 chip barker code to demonstrate my first statement. This Barker code has a similar property to C/A codes in that a correlation with shifted versions with itself (autocorrelation) will have a very low correlation relative to when the code is aligned.

1 0 1 1 0 1 1 1 0 0 0

This code written with 1's and 0's has a mean of approximately 1/2, so you can optionally remove the mean and similarly for convenience scale by a factor of two by mapping 0 to -1 such that the code is:

1 -1 1 1 -1 1 1 1 -1 -1 1

This results in the following correlation when the sequences are aligned in time (which could be used in direct-sequence-spread-spectrum, DSSS, to send a data symbol with value = "1"):

And the following if the code was inverted (which could be used in DSSS to send a data symbol with value = "0"):

And the following for any other shift (I did a rotational shift of one sample here, corresponding to m = -1):

With the following autocorrelation result as a plot:

We could similarly instead perform the correlation with the following operations for multiplication (basically XOR function, could do the same with XNOR and get an inverterted result):

1 x 0 = 0

0 x 1 = 1

1 x 1 = 1

0 x 0 = 0

Resulting in a similar result with an offset and scaling in the correlation output instead of an offset and scaling the correlation input. So in this case the maximum correlation would be 11, the maximum correlation when inverted is 0 and the autocorrelation for any other offset would be 5. Double this and subtract 11 and you get the same result as before.

(Note unlike maximum-length pseudo-random sequences which will also cross-correlate to -1 at all offsets when you use +1 and -1 values, C/A codes have cross-correlation values of either 1, +63 or -65, which is still significantly less than the 1023 value when aligned).

The reason this works for GPS in the case of low SNR is when you add each of the 1023 chips with noise present, if the noise on each chip is independent ("white noise"), then the signal magnitude will accumulate by a factor of $1023$, but the standard deviation of the noise will only accumulate by a factor of $\sqrt{1023}$ thus in terms of SNR you get a $10Log_{10}(N)$ processing gain. Consider a simpler case of adding two independent Gaussian random variables with equal distribution and equal but non-zero mean: The mean will double but the standard deviation will only increase by the $\sqrt{2}$.


Computation of the Inverse DCT (IDCT) Using DCT Or IFFT Is there a way to compute the inverse discrete cosine transform (type-2) by leveraging either a DCT, FFT, or IFFT algorithm? I have seen ways to compute the DCT using FFTs, and I've seen ways to compute IFFT using FFT. I can't quite find a simple example with description.

Have a look at Fast DCT Algorithm (PDF Version).

It has both DCT and Inverse DCT using DFT (FFT).
They show how to do a DCT and Inverse without the reflection trick.

The standard (Less efficient then above) way doing so is:

numElements = 10;

vX = randn(1, numElements);

disp(vX); %<! Diusplay the input signal

vDctRef = dct(vX);

% Forward DCR using FFT
vXX     = [fliplr(vX), vX]; %<! Mirroring
vXDft   = fft(vXX);

vGrid = [0:((2 * numElements) - 1)];

vShiftGrid = exp(-1j * 2 * pi * (numElements - 0.5) * vGrid / (2 * numElements));

vXDft2 = real(vXDft ./ vShiftGrid);

vDct    = vXDft2(1:numElements) / sqrt(2 * numElements);
vDct(1) = vDct(1) / sqrt(2);

disp(vDctRef)
disp(vDct)

% Inverse DCT Using FFT
vDct(1) = vDct(1) * sqrt(2);
vDct    = vDct * sqrt(2 * numElements);

vXDft = [vDct, 0, -fliplr(vDct(2:numElements))] .* vShiftGrid;
vXX = ifft(vXDft);

vX = real(fliplr(vXX(1:numElements)));
disp(vX);

A reference code for the IDCT by the more efficient method (As in reference):

numElements = 10;

% Signal (DCT Coefficients)
vXDct = randn(1, numElements);

% Reference Inverse IDCT
vXRef = idct(vXDct);

% Inverse IDCT Using FFT
vGrid = [0:(numElements - 1)];

vShiftGrid      = exp((1j * pi * vGrid) / (2 * numElements));
vShiftGrid      = vShiftGrid * sqrt(2 * numElements);
vShiftGrid(1)   = vShiftGrid(1) / sqrt(2);

vTmp = vShiftGrid .* vXDct ;
vTmp = real(ifft(vTmp));

vX = zeros(1, numElements);

for ii = 0:((numElements / 2) - 1)
    vX((2 * ii) + 1) = vTmp(ii + 1) ;
    vX((2 * ii) + 2) = vTmp(numElements - ii) ;
end

disp(vXRef);
disp(vX);

Another reference is given by Lecture Notes for the Course MAT-INF2360 - Fourier Theory, Wavelet Analysis and Non Linear Optimization.
Have a look at 4.2 Efficient implementations of the DCT at page 115.

Enjoy...

Remark
When it is written DCT above it refers to DCT Type II.


Sampled AC Voltage measurement I want to implement a AC Voltage measurement function for a DMM. That DMM can sample up to 1,000,000 Samples/s.
With the sampling function I wanted to sample an input sigal (e.g. sine signal) and measure/calculate the RMS value.
Since the formula for the RMS calculation in known, that part is not a problem.

My question is, how do I choose the right sampling rate and the number of samples? For different signal frequencies I need different sampling rates and different sample sizes (considering  the Shannon Theorem and the duration of the sampling process). In my case, I'd like to be able to measure the RMS of signal from 50Hz up to 500kHz.

Is there a certain standard how to implement such a function?

RMS is essentially a three step process on a signal:  Square it, mean it, root it.

Squaring is not so difficult in software.  Two things to keep in mind here.  First, squaring fixed point data will double the data width, so make sure you have enough precision to avoid overflow.  Second, squaring a signal will double its bandwidth.  Make sure you have enough margin in your sample rate to either eliminate aliasing, or that the aliasing doesn't affect your baseband signal.

A square root can be computationally expensive, but it isn't exactly difficult to do.  I does have an impact on bandwidth, but this could be treated as a data representation problem, for which you would not consider it to affect bandwidth.  I wouldn't worry about it.

A mean technically requires the full data set, so it doesn't apply to real time processing.  To get around this, we can instead apply low pass filtering.  Low pass filtering is a big enough topic that I don't think it warrants much discussion here, but you could any standard filtering techniques, such as IIR, FIR.  I would recommend trying some and maybe asking another question on here if you need more information.

It isn't explicitly required, but it is likely that you would need to decimate your signal at some point.  For example, your sample rate may be 1MHz, but the user interface might only refresh at 100Hz.  This would happen either after or during your filtering stage.  You could apply a pair or second order maximally flat IIR low pass filters with a cutoff of say 25Hz, then throw away 9,999 out of every 10,000 samples.  You could also save 10,000 samples in a buffer at a time, and calculate the average of this set.  This is (I believe) how an oscilloscope does its math, but keep in mind that the accuracy wouldn't be great, because the filter is not aggressive enough for the decimation.  You could 'smooth' the data by applying a window to the buffered data, but you would still get some aliasing if your buffers do not overlap.

Lastly, the bandwidth of the input signal.  At 1MHz, 500kHz is at Nyquist, so it is the absolute maximum bandwidth you can sample without aliasing.  However, this doesn't provide any margin, so your sampling is likely to be sub par because of real world sampling limitations, such as quantization error, and noise.  Additionally, squaring a 500kHz would create content at 0Hz as well as 1MHz.  The 1MHz content would alias back down to 0Hz, causing distortion of the 0Hz signal.  As such, you would want to build in some margin to avoid aliasing.

All this being said, I don't know how much proper signal processing technique goes into your off the shelf DMM.  In fact, it is possible to do a fair bit of this process in the analog domain to offload processing power.  It seems more likely that someone designing a DMM would be inclined to throw high precision, high bandwidth ADCs at the problem and brute force a solution that provides a number.  That number might be accurate for a steady state sine wave at a given frequency range, but inaccurate for real world signals.  However, it is unlikely that the user would either notice or care.


Harmonics in Quantal Histogram I'm studying what I believe to be a system that produces a (noisy) quantal response (i.e. there should be a minimum increment of response). To do this, I've plotted a histogram of some ~6000 amplitude values from the system into small bins, and performed an FFT on the histogram, expecting that a common quantal value should appear as the lowest frequency "spike". However, my spectral response:
 looks odd to me. What I see are what look like harmonics.

Why does this occur in the quantal release problem? Is this lowest frequency actually my quantal response, or is it a typical electrical response/artefact from the FFT of a signal with so many discontinuities? I thought that odd-numbered harmonics are the most powerful, but here they seem to be linearly declining. Why am I seeing other smaller spikes appear 1/3rd of the way through each of the "main" harmonic bands?

The key problem with this histogram is that it looks more like a collection of impulses rather than a "waveform" and this "shows" at the output of the Discrete Fourier Transform (DFT).

If you have to use the DFT then you need to preprocess the histogram so that it becomes smoother. You can do this with a moving average filter that will effectively integrate the impulses into smooth "humps". The downside of this is that you will have to tune the integrating period of the moving average filter to something that makes sense for your problem. And while doing this you will also notice that the waveform might change a bit from the expected depending on the time scale (or in the case of the histogram the resolution) that results from the use of a filter.

You can then pass the smoothed waveform to the FFT or, if your objective is to determine modes occuring at regular intervals, pass the smoothed histogram "waveform" to something like Kernel Density Estimation which seems more suited to your application.

Kernel Density Estimation will try to approximate the smoothed histogram output as a sum of gaussians with different widths and centers. Ideally, if these modes resemble gaussians, you might get one gaussian per mode you are trying to locate and then the regularity of the modes will show up in the centers of the fitted gaussians.

Hope this helps


State space system identification Suppose I have a real (physical) dynamical system with some sensors and actuators, and I also have an idealized state-space model of this same system.  How, in general, can I adjust the model to match the real system?

My experience is with frequency-domain analysis.  My first reaction would be to measure swept-sine transfer functions, and to calculate parameters of my system from the location of resonances or poles in the measured data.

Is there a better, more general technique for system identification in the state-space representation?

The experimental frequency domain data gives you information on the input-output behavior of the system. If the system is assumed to be linear and time-invariant (LTI), this information will be in the form of a transfer function. Given a transfer function, there are infinitely many state-space realizations that yield the same input-output behavior and, thus, are I/O-equivalent.

Having said that, you may want to take a look at the Adaptive Control literature.


why bandpass filtering (butterworth) centers the signal I run bandpass filtering via the butter() in matlab to apply butterworth bandpass filtering to find out the signals are centered at zero after the filter. Would this always happen for any kind of low/high-pass filtering and why?

High pass filtering means that low-frequencies, including DC, are removed. So if your signal has a DC offset (is not zero mean), then this will be removed by a highpass filter or a bandpass filter that does not include 0 in the passband.


To implement successive interference cancellation (SIC) process, should the receiver know the power allocated to each symbol? My question is as follows:

To implement successive interference cancellation (SIC) process in downlink non-orthogonal multiple access (NOMA), should receivers know all power allocated to each symbol?

Consider a single-cell network consisting of one base station (BS) and $N$ users.
Let $x_i$ be the symbol for user $i$, $h_i$ be the complex-valued channel gain between user $i$ and the BS, $p_i$ be the power allocated to symbol $x_i$.
Then, the transmitted signal is given by
$$
x = \sum_{i=1}^N \sqrt{p_i} x_i.
$$
The corresponding received signal at user $k$ is given by
$$
y = h_k \sum_{i=1}^N \sqrt{p_i} x_i + n_k,
$$
where $n_k\sim\mathcal{CN}(0,\sigma^2)$ is the additive white Gaussian noise of user $k$.
Assuming that $\lvert h_1 \rvert^2 \le \lvert h_2 \rvert^2 \le \cdots \le \lvert h_N \rvert^2$, user $k$ decodes its own symbol $x_k$ according to the following process:

Decode $x_1$ by treating $\sum_{i=2}^N h_k\sqrt{p_i}x_i + n_k$ as noise.
Subtract $\hat{x}_1 \triangleq h_k\sqrt{p_1}x_1$ from $y$.
Decode $x_2$ by treating $\sum_{i=3}^N h_k\sqrt{p_i}x_i + n_k$ as noise.
Subtract $\hat{x}_2 \triangleq h_k\sqrt{p_2}x_2$ from $y-\hat{x}_1$.
Repeat to (i) decode $x_l$ by treating $\sum_{i=l+1}^N h_k\sqrt{p_i}x_i+n_k$ as noise, and (ii) subtract $\hat{x_l} \triangleq h_k\sqrt{p_l}x_l$ from $y-\sum_{i=1}^{l-1}\hat{x}_i$ until $x_k$ is decoded.

In my opinion, to decode its own symbol $x_k$, user $k$ should know $p_l$ for any other user $l$ such that $\lvert h_l \rvert^2 < \lvert h_k \rvert^2$.
Is this opinion true? If this is true, the BS should inform the user about the power allocated to each symbol. Nevertheless, is NOMA useful?

Yes, one would require to know the power of other users. There are two options, the network can either signal the UEs this information in a control frame or the UEs need to estimate the power of the streams allocated for other users (knowing its own power allocated and some smart signal processing).
In any case interference cancellation techniques are suspetible to the "effectiveness" of the cancellation. Usually the SNR degardes rapdily with decrease in estimation accuracy or error in knowledge of the parameters of the other users, such as power


How to derive low pass filter $\frac{N}{N+S}$ I watched the video for PID control system where it mentioned the Laplace domain function and low pass filter to be $\frac{N}{N+S}$.
I used asymptotic analysis to see that it made sense. However, I'm interested in how it was originally derived. How did Laplace domain function lead to the low pass filter?

I assume that you want to see that the transfer-function magnitude is inverse proportional to the frequency. In this case, you replace $S$ with $i\omega$ and you get that the transfer function magnitude decays with the frequency$|G(\omega)|=\left|\frac{n}{ n+i\omega}\right|$


Extract repeating structure from an image I want to find repeating structures from the image .But I have no idea how to do it.Please help me.Actually i want to find how many lines and oval shapes look like similar and mark these by separate color.

You want to detect two things: lines/curves and circles. There are well-known algorithms to do both, so let's find them separately. I will present a solution in Mathematica. You can translate to your language.

The first step is to eliminate the background. We do this by finding the background:
img = Import["http://i.stack.imgur.com/Olczj.jpg"]
bimg = ColorConvert[img, "Grayscale"]
Blur[Closing[bimg, 20], 15]

Then dividing by it:

  evenImg = Image[ImageData@bimg/ImageData@Blur[Closing[bimg, 20], 10]]

Detect the edges:

  binImg = EdgeDetect[evenImg, 1, .1]

Now we just need to fill in the strokes:

  smImg = ImageResize[Dilation[binImg, 2], {500}]

We are finally ready for some detection. Let's start with the lines/curves:

  Show[smImg,
   Graphics[{Yellow, Thick, Line /@ ImageLines[smImg, 0.25]}]]

What about the circles? Let's start by looking at all the connected components:

  morph = MorphologicalComponents@smImg; Colorize@morph

The circles are there, all right. We just need to filter out the objects that are not circular. We can do this by looking at the circularity, elongation, eccentricity, or Euler characteristic. Circles are easy to distinguish from other shapes by these properties.

noncircles = #[[1]] & /@ ComponentMeasurements[morph, "Elongation", # > 0.3 &]
circles = morph /. x_ :> 0 /; MemberQ[noncircles, x] //
  Colorize // RemoveBackground


About Fourier transform of periodic signal In Fourier transform for periodic signal, I checked different books and I found a different explanation in each book. Let's take the explanation in Signals and Systems by Rajeshwari & Rao:

  The resulting Fourier transform for a periodic signal consist of a train of impulses in frequency, with areas of impulses proportional to the Fourier series coefficients.

  To suggest the general result, let us consider $x(t)$ with Fourier transform $X(\omega)$ which is a single impulse of area $2\pi$ at $\omega=\omega_0$, that is,  $$X(\omega)=2\pi\delta(\omega-\omega_0)$$
  To determine $x(t)$ for which this is Fourier transform we can apply  the inverse Fourier transform to obtain
  $$ x(t) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} 2\pi \delta(\omega-\omega_0)e^{j\omega t}d\omega$$

The things I want to ask is:

If we have Fourier series of a periodic signal which will be one impulse, then the Fourier transform of that impulse will be the same single impulse? That's what it is explained above?
Why we used shifted impulse? Why we can't take $\delta(\omega)$

We can figure out what's going on if we first understand a simple identity and then just compute the Fourier transform of the periodic function.

A useful identity

First let's prove that
$$D(\omega - \omega') \equiv \int_{-\infty}^\infty dt \, e^{i (\omega-\omega') t} = 2\pi \, \delta(\omega - \omega')\,. $$
We just use a test function $\tilde{g}(\omega)$:
\begin{align}
\int_{-\infty}^\infty \frac{d\omega}{2\pi}
\underbrace{ \left(\int_{-\infty}^\infty dt \, e^{i (\omega-\omega') t} \right)}_{D(\omega - \omega')} \tilde{g}(\omega)
&= \int_{-\infty}^\infty dt \, e^{-i \omega' t}\left( \int_{-\infty}^\infty \frac{d\omega}{2\pi} \tilde{g}(\omega) e^{i \omega t}\right) \\
&= \int_{-\infty}^\infty dt \, g(t) e^{-i \omega' t} \\
&\equiv \tilde{g}(\omega') \, . \\
\end{align}
This shows that by definition $D(\omega-\omega') = 2\pi \, \delta(\omega - \omega')$.
We didn't have to say anything weird about distributions or use any magic.

Fourier transform of periodic function

Now consider a periodic function $x(t)$ with period $T$.
Since $x$ is periodic we can write it as a Fourier series
$$x(t) = \sum_{n=-\infty}^\infty \hat{x}_n e^{i 2\pi n t/T} \, .$$

Now let's compute the Fourier transform,
\begin{align}
\tilde{x}(\omega)
&\equiv \int_{-\infty}^\infty x(t)e^{-i\omega t} \, dt\\
\text{plug in Fourier series} \qquad &= \sum_{n=-\infty}^\infty \hat{x}_n \underbrace{\int_{-\infty}^\infty dt \, e^{i(2\pi n/T -\omega) t}}_{D(2\pi n / T - \omega)} \\
&= \sum_{n=-\infty}^\infty 2\pi \, \hat{x}_n \delta(\omega - 2\pi n / T) \, .
\end{align}
So there we go: the Fourier transform of a periodic function $x(t)$ is a sum of delta functions appearing at $x(t)$'s Fourier series frequencies (i.e. the frequencies $2 \pi n / T$).
The weight of each delta function is $2\pi$ times the Fourier series coefficient $\hat{x}_n$.
The factor of $2\pi$ is there simply because we used the Fourier transform convention with angular frequency $\omega$ but the Fourier series convention with regular frequency.


Generating QPSK Constellation Diagram I was wondering how a diagram like this is generated.

I am generating a QPSK signal using MATLAB by taking a cos and sin wave, multiplying each one by either 1 or -1 and adding them together to form a real valued QPSK signal. I then simulate an AWGN channel and send it to my demodulator. I mix it down using another cos and sin wave plus Low Pass Filter to get my recovered I and Q data. At this point I am left with the baseband which I can recover my data from by taking the angle of the I and Q components.

I am wondering at what point in that chain I could generate a the plot I post above.

Thank you

This is what I do when making such a plot for a QPSK modulation:

First I upsample the I and Q data so that I can see more of the trajectory between samples given (assuming you start with at least 2 samples per symbol otherwise, due to Nyquist, the information on the trajectory between symbols is not available, but as long as you have 2 or more samples, all the information of what happens in the first Nyquist zone is preserved):

Is = resample(I, 10, 1);

Qs = resample(Q, 10, 1);

and then I make a plot of I vs Q

plot(Is,Qs)

I used a resampling increase of 10, but you could make it higher if desired. The result is the interpolated samples for the full trajectory during symbol transitions.


Representation, Bandwidth and composition of Digital Signals I have started learning data communication and I'm not able to grasp certain topics.

First

(i) I know that analog signals are continuous and digital signals are discrete. Analog signal are represented by a continuous time-varying quantity like voltage. How are digital signals represented in actual practice? I mean to ask how are the 0s & 1s actually represented/transmitted over a transmission medium? Are they represented by a sudden spike then staying stable and then a sudden fall of the time-varying quantity like voltage?

Second,

(i) I know bandwidth is the range of frequencies that make up a composite analog signal. So, if there are two analog signal with frequency 5 Hz and 10 Hz making up a composite signal, will it's bandwidth be (10-5=) 5 hz?

(ii) When we say bandwidth of a transmission medium is 5 hertz, do we mean to say that it allows any composite signal to pass through it, in which difference between the upper frequency and lower frequency making up the composite signal is 5 Hz?

(iii) Can we define bandwidth of a single sine analog signal?

Third

(i) Are digital signal nothing but just superimposition of theoretically an infinite number of analog signals? If yes, then why and how?

(ii) How do we define frequency and bandwidth(in Hertz, not bps) of a digital signal?

Sorry if these questions seem a bit too easy for many but as a self-learner with no college help, It's quite tough for me.
Thanks in advance.

First

  (i) I know that analog signals are continuous and digital signals are
  discrete. Analog signal are represented by a continuous time-varying
  quantity like voltage. How are digital signals represented in actual
  practice? I mean to ask how are the 0s & 1s actually
  represented/transmitted over a transmission medium? Are they
  represented by a sudden spike then staying stable and then a sudden
  fall of the time-varying quantity like voltage?

Physical representation of discrete-time (or digital) signals is also through analog waveforms, but the information they carry is discrete. They can be transmitted like pulses of $5$ V, $0 $ V in a digital computer, or sine waves of different frequencies like in a digital telephone line.

  Second,

  (i) I know bandwidth is the range of frequencies that make up a
  composite analog signal. So, if there are two analog signal with
  frequency 5 Hz and 10 Hz making up a composite signal, will it's
  bandwidth be (10-5=) 5 hz?

If you include all frequencies from $0$ Hz upto the max in the signal, then the bandwidth of your composite baseband signal will be $10-0 = 10$ Hz., counting only the positive excursion of frequencies. If you want to consider your signal as bandpass then its bandwidth will be $10-5=5$ Hz.

  (ii) When we say bandwidth of a transmission medium is 5 hertz, do we
  mean to say that it allows any composite signal to pass through it, in
  which difference between the upper frequency and lower frequency
  making up the composite signal is 5 Hz?

No. The actual frequency range of that passband is also important; a medium that can transmit a band in [15-20] Hz cannot transmit a signal in the band [85-90] Hz, even though both have the same bandwidth. You should shift the signal spectrum into the channel band for transmission.

  (iii) Can we define bandwidth of a single sine analog signal?

Bandwidth of $x(t)=\cos(2 \pi f_0 t)$ is its frequency $f_0$ if you treat it like a baseband signal, but otherwise $0$ Hz if you treat it like bandpass signal...

  Third

  (i) Are digital signal nothing but just superimposition of
  theoretically an infinite number of analog signals? If yes, then why
  and how?

digital signal is a misnomer here. It's true that an ideal pulse waveform requires infinite spectrum and infinite number of sine waves to construct it, but a pulse is not the only way to carry digital information as the answer to the first question already indicated...

  (ii) How do we define frequency and bandwidth(in Hertz, not bps) of a
  digital signal?

digital signal is a misnomer here, it's the discrete information that you actually should care, that information is quantified by an entropy measure. Then this will be transmitted via line signalling with analog bandwidth in Hz depending on the efficiency of signalling waveforms.


A question on Fourier Series and the frequency of the sinusoids On studying about Fourier series, I encountered 2 doubts:

How is it that a non-periodic function has a Fourier series?
When expressing a periodic function as summation of sinusoids, why is the fundamental (minimum) frequency taken to be that of the periodic signal? Why cannot one of the sinusoids have a frequency less than the fundamental?

I mean, we think of representing the signal in terms of sinusoids with frequencies greater than and equal to the fundamental frequency, and of course, there is the DC part as well. But, why is there no sinusoid with frequency less than the fundamental as well?

Is there any mathematical logic for the second point?

A function that is not periodic doesn't have a Fourier series. What you can do is take a finite length interval and construct a Fourier series which approximates the function in that interval (and continues it periodically outside that interval).
The sum of several sinusoids with one of the sinusoids having (fundamental) period $T$ can never have a period less than $T$, simply because that one sinusoid will destroy periodicity (with period $<T$) due to its larger period.


Difference between ycbcr and ycbycr I see in a paper mentioning ycbcr as well as ycbycr. I am wondering what is the difference between the two? In the figure, the first two lines show ycbycr but the two bottom lines show ycbycr?

Brightness (Y) component is sampled at a twice rate compared to the color components in YCbYCr coding. Such coding makes use of a fact that human eye is more sensitive to brightness than to color and thus helps to limit bitrate while maintaining acceptable visual quality. Decimated color component samples are reconstructed by interpolation. Alternatively, YCbCr honestly transmits color components for each pixel.

That can be compared to a well-known Bayer mask that samples G component twice as fast compared to R and B because human eye is more sensitive to it.


Improving the communication speed on a 10km cable (with high capacitance) I'm working on improving the modem communication speed that can be achieved over a 10km cable. The current modem achieves a speed of 9600 bps via hardware (phase locked loop) and I'm looking for ways to at least double this, and do it all in firmware. Here are some specs for the project.

Cable characteristics

Length: 10,000m
Capacitance: 1.48μF (148pF/M)
Resistance: 320 Ω (32 Ω/km)

Driver signal

10 Vpp (+10 V to -10 V) @ 250 mA
Can be used for Amplitude and Frequency modulation
Currently cannot perform phase modulation

Available bandwidth on cable

30 KHz to 100 KHz
Attenuation: 100KHz signal results in a third of the amplitude of the 30KHz signal

Current method of modem communications (want to improve)

A hardware implemented Phase Locked Loop (PLL)
Achieves a speed of 9600 bps

The aim is to improve the 9600 bps speed and also to implement the modem in firmware. The firmware target is an ARM based 32-bit microcontroller - Infineon's XMC4500 with onboard ADC (firmware to be written in C). I have a working prototype, but the best I can achieve so far is 10,500 bps which is only a small improvement over the current system. The method used was to implement FSK demodulation using the Goertzel algorithm. I was detecting 4 frequencies in the range of 40KHz to 80Khz, which gives 2 bits per signalling period. On the microcontroller performance, I've found that 4 frequencies is the maximum number that can be attempted using the Goertzel method. Also the microcontroller ADC's practical maximum sampling rate is 700 KHz.

Does anyone have a suggestion how to achieve a decent speed improvement using frequency and/or amplitude modulation with such limitations as these? I'm new to DSP so I've likely missed some other better techniques that could be used.

Many thanks in advance.

PROGRESS UPDATE (10/10/2013)

Since posting this question and after some simulation and testing, the communications scheme has been changed to be asymmetrical since only the downlink side actually needed to be upgraded. The uplink speed remains as is - a PLL running at 9600 bps, but the downlink speed ends up being a faster ~64,000 bps, with the downlink demodulator hardware running on something equivalent to PC level.

This hardware upgrade enables the use of an FFT and an 8-bit symbol using 8 frequencies. The FFT parameters being used are a sample rate of 512,000 SPS and a number of bins of 64. The frequencies being detected are in the range of 32,000 Hz and upwards, with a spacing of 8,000 Hz. The ADC is 12-bit. I've also implemented Reed-Solomon forward error correction, which may be used to reduce error rates.

It seems your hardware is limited. If you can't detect more finite frequencies and you can't implement another method of modulation you might have reached your maximum in physical bits with your hardware.

May I suggest a bit of coding theory? Implement LZW compression. The data rate will vary and it might not give you double the speed but will definitely increase the data rate.

Its already produced in c, try it out.

http://marknelson.us/attachments/1989/lzw-data-compression/lzw.c


Peak Detection Of Every Pulse I have this signal as the following figure:

Now, I need to know each beginning index of each pulse as described in the following figure:

I have tried to detect the peaks through obtaining the cross correlation function with the sequence that I am sending and the Signal that I am receiving (which is exactly as if am doing auto-correlation). But they are not aligned as the following figure:

Can you suggest a solution to this problem, please? I am kinda stuck and this is my last step in my project.

I can show the codes if needed.

filtering -> thresholding -> edge detection

apply low pass filter:

averageFilterCoefficient = 0.01
for each sample of x:
  y[i] = y[i-1]*(1-averageFilterCoefficient) + x[i]*averageFilterCoefficient

this simple code should remove all high frequencies and produce more smooth signal. Now apply thresholding:

threshold = 100
for each sample of y:
  if y[i] > threshold then
    z[i] = 1
  else
    z[i] = 0

Now you have nice logical signal, lets find rising edges:

for each sample of z:
  if z[i-1] = 0 and z[i] = 1 then
    print i

You just need to find the right values for exponential average filter
coefficient and thresholding level.


Fourier transform as the integral of a parameter multiplied by an homogeneous wave Can a Fourier transform in space be interpreted as the integral of a parameter multiplied by an homogeneous wave $\sigma$?
where $\sigma$ is:
$\sigma$=$e^{-ikx}$
Are there papers or book that illustrates this interpretation?

given your definition, and the definition of the Fourier transfer:
yes. That is literally the definition:
$$\int_{\mathbb R} s(x) e^{i2\pi f x} \mathrm dx$$
By comparison, with $k= 2\pi f$, you get the Fourier transform. Every book uses that definition, so you'll find this "interpretation" in every book that introduces the Fourier transform. (This really feels like exactly the same $k$-space / location transform I learned in my first solid state physics lecture, so I assume every modern solid state course for physicists does pretty much that, but more in-depth.)


Calculating the Spectral Centroid of a Signal Here is the equation:

$$C_i=\frac{\displaystyle \sum_{k=1}^{W_{f_L}}kX_i(k)}{\displaystyle\sum_{k=1}^{W_{f_L}}X_i(k)} $$

The MATLAB code for the equation is:

sum(m.*window) / sum(window);

Let's say length(window) = 882 that contains the first half (abs) of the magnitude spectrum, fs = 44100.

Now calculate m:

m = ( (fs/(2*window)) * [1:window] )';

The length(m) = 882 and increments in 25, which is the sample range (e.g. [25 50 75 100 ... 22050]).

I'm trying to determine a better insight in terms of what is happening when both terms are being divided.  It is clear that each element of m is multiplying each element of window, but what is actually happening?  The book I'm reading just states that the Spectral Centroid is the center of 'gravity' of the spectrum.

The function you wrote is basically like calculating the center of mass of a body or for that matter the Mean Value of a Distribution.

If you remember from Probability, the Mean is given by:

$$ E \left[ y \right] = \sum_i {y}_{i} p \left( {y}_{i} \right) $$

Now, let's think of the Frequency Function (The DFT) as a probability function.
In order to do so we'll normalize it to have a sum of 1.

This is done by dividing each element by the sum of elements.
Then just apply the Mean operator to calculate the Mean of this "Distribution":

$$ {C}_{i} = \sum_{k} \underset{{y}_{i}}{\underbrace{k}} \underset{ p \left( {y}_{i} \right ) }{\underbrace{ \frac{{X}_{i} \left( k \right)}{\sum_{k} {X}_{i}
 \left( k \right) } }} $$

The result you'll get is exactly what you have written above.


Advantages of 2-VSB vs BPSK? Both BPSK and 2-VSB transmit one bit of information per smybol and both modulate the real part of the signal.
The difference is, as far as I understand, that BPSK does not modulate the imaginary part of the complex signal at all, while 2-VSB uses hilbert transformation to suppress the lower sideband of the transmitted signal.

Are there any advantages of doing this or is it just adding complexity without any advantages?

Suppressing the lower sideband decreases the bandwidth by a factor of $2$ compared to BPSK. Since we're not talking about single-sideband (SSB) modulation, but vestigial sideband (VSB) modulation, we don't gain a factor of $2$, but a factor that is slightly smaller than two, because a part of one of the sidebands is retained.

Note that both BPSK and 2-VSB carry the information only in the (bi-polar) amplitude of the signal, i.e., only in the in-phase component, not in the quadrature component. The complex baseband signal for 2-VSB only occurs because the larger part of one of the sidebands is filtered out. This, however, is not necessarily done by complex baseband processing (using the Hilbert transform), but it can be done by applying a band pass filter after modulating the baseband signal. The latter option only uses real-valued processing.


GNU radio, FSK channel estimator A quick question, is there any way to do channel estimation on a (G)FSK receiver. I'm trying to build a DASH7 receiver.

I found that there is a block for a OFDM channel estimator, but i can't find it for FSK.

In FSK, you don't care about the frequency-varying phase properties of your channel – they simply don't matter to your signal, so a (relative) amplitude estimate for all mark frequencies would technically suffice.

Most FSK receivers don't even do that – it's only helpful to know that "mark frequency A gets through stronger than frequency B" if you use that to adjust your decision boundaries (e.g. for an MAP decision, you'd need that knowledge).

I don't have access to any DASH7 phy layer documentation, but if it's really derived from RFID tag technology, then you can be 100% certain the thing works without channel state information. Complexity wouldn't allow that. That doesn't mean you can't make it work better at further distances – but it does imply that if you're well within the specified ranges and environments, getting a channel estimate will not help you significantly much.
Wikipedia says "range up to 2 km, less than 200 kHz bandwidth", and that basically means your second observable multipath component would have to have a 1.5 km longer way to travel than the first one you receive – that is physically practically irrelevant, as that delayed version of the signal will be so much weaker than the quicker path that it can't lead to significant frequency-selectivity.

So the answer to your question is "you don't do that at all".

The rest of this answer, hence, is of merely academic interest:

Now, if your channel is actually strongly frequency-selective, so that knowing the relative amplitude becomes desirable, you can also be sure that the transition between different frequencies (which is Gauss shaped, transmitter-side) isn't nice anymore – and that's bad for your Gauss-expecting receiver.

Now, the only thing that could be FSK-specifically done here is some form of nonlinear emphasis on some values – that really doesn't sound very tempting, because estimating the necessary distortion to counteract the nonlinear memory-affected effect of a frequency selective channel on the FM demodulation of a frequency transition sounds pretty hard. And whenever possible, we try to compensate systems where they are linear – that gives us at least some guarantees about invertibility/well-posedness of the inverse problem.

So, don't try to deal with your FSK in any special way. Instead, just do what works for any transmission type:

transmit something you know – for example, a preamble, or a known data packet – that is as white as you get with the given scheme (in this case: transmit pseudorandom white data through your GFSK transmitter)
Correlate with that at the receiver
Since that correlation will be mathematically equivalent to (noisy) convolution of your transmit signal's autocorrelation function and the channel impulse response: It's your channel estimate. You don't need it any better than that, because it excites exactly the frequencies present in your signal.

The first step – transmission of something you know – is in many systems done via a decision feedback equalizer, where "something you know" is what you received and decided for, but that means block-wise operation in an FSK case, since your system is not linear in modulation pulses, but has memory.

TL;DR: Don't. Your system design implies the channel is flat, and hence you don't need a channel estimate. Spend your time on something else!


IQ data QPSK Demodulation I have an IQ data of QPSK signal , when I apply a PLL loop I obtain constellation in figure below. What is this phenomenon, and how can I obtain the good constellation??

The second figure is about variation angles, is there any idea about the red countoured effect?

after adding symbol syncronisation , i obtain constellation below , but there  is also sommes error that causes incorrect demapping ;; is there any suggestion ????

I suspect you're not doing symbol synchronization. You need to synchronize and then sample the output of your receiver filter at the correct times. Only then will you see the tradiational constellation diagram. Look for "Gardner algorithm" on this website; another well-known algorithm is Mueller and Muller.

One thing you can do to get confidence that your signal is correct is to plot its eye diagram. If the eye is open, then all you need to do is synchronize and sample. If it's closed, you have a different problem.


Getting Composition of a Mixture of Colors I wanted to validate a thought on blending of colours digitally.

Imagine I have two colours, $C_1$ $(R_1,G_1,B_1)$ and $C_2$ $(R_2,G_2,B_2)$. Now, I create colour $C3$ as follows:
$C_3 = \alpha C_1 + (1-\alpha) C_2$

Where $\alpha$ is a value between $0$ and $1$.
So, yes, basically, I am trying to put one object over another with a transparency associated with it.

What I want to understand is that, if I have $C_3$, $\alpha$ and $C_1$ in the above equation, will I get an accurate value of $C_2$?

Math suggests that I should, but I just want to confirm that I am not missing anything.
Also, can I keep a really high but $< 1$ value for $\alpha$ and still accurately get $C_2$?

You actually have 3 independent equations.
If one of them is solvable, all three are.

Now, the equation is given by:

$$ z = \alpha x + \left( 1 - \alpha \right) y $$

Assume $ \alpha, x, z $ are known, then:

$$ y = \frac{z - \alpha x}{1 - \alpha} $$

For any value of $ \alpha \neq 1 $ the solution is valid.
Moreover, if $ alpha = 1 $ then $ z = x $ and $ y $ is irrelevant.

Looking at the other case, restoring $ x $, yields:

$$ x = \frac{z - \left( 1- \alpha \right) y}{\alpha} $$

For any value of $ \alpha \neq 0 $ the solution is valid.
Moreover, if $ alpha = 0 $ then $ z = y $ and $ x $ is irrelevant.

As you wrote, in order to have a Linear Interpolation the parameter $ \alpha $ must be in the range $ \left[ 0, 1 \right] $.


time domain data into FFT using Matlab I do like to generate the FFT from time domain history data (experimental lab) using MATLAB. I am not so sure if my coding is right? Can someone verify the coding.

My code is as follow:

 clc;
 signalx = Ch1;
 N=6678;
 t=Time;

 %Plot FFT Vehicle signal
 zx = fft(signalx);
 plot(abs(zx));
 title('FFT Vehicle');
 xlabel('Frequency [Hertz]')
 ylabel('Magnitude');

Just to expand what I wrote in the comment. This is how you should work with the FFT in Matlab.

Compute the FFT of your time signal:

NFFT = 2^nextpow2(N);
Y = fft(y,NFFT)/N; % the division by N is to scale the amplitude

Then compute the sampling frequency: (EDIT skip this 2 lines of code if you already know the value of $Fs$)

dt = t(2)-t(1);
Fs = 1/dt;

This works if your data are sampled at constant sampling frequency. If it is not the case you need to resample them (it is very easy to do if you convert your data into timeseries and then use the resample function for timeseries. Look here for more details: What is an algorithm to re-sample from a variable rate to a fixed rate?)

Finally, plot the result:

f_plot = Fs/2*linspace(0,1,NFFT/2+1); % vector of frequencies for the plot
figure()
plot(f_plot,2*abs(Y(1:NFFT/2+1))) % multiplicated by 2 to recover the energy related
% to the negative frequencies since in this way only half spectrum is plotted
xlabel(’Frequency’)


Acoustic record and store devices I am working on a project involving analysis of bird sounds and I want to record these sounds in slightly remote locations. It is enough if the device can only record and store them (maybe in an SD card). Are there any devices available that I can readily use or what would be the best way to do this ?

There are plenty of devices that can do this. Simply Google portable recorders. They generally all have a microphone built in. Some have line level inputs, and some can also take XLR inputs and provide phantom power for condenser microphones.

If you don't need the XLR inputs, I can recommend the Zoom H1. It has a stereo microphone, can also take a line level input and is fairly cheap. I have used this to record various foley sounds including bird songs. The sounds save onto an SD card in high quality 24-bit/96 kHz WAV files for analysis on a computer later.

It's best to get a little tripod for it too to stop any handling noise.


How to calculate the frequency of a signal without knowing the sampling frequency using Matlab I have a signal and I am using Matlab command pwelch to calculate the frequency of the signal, but the frequency I obtained is changed as I change the sampling frequency.

For example, when using sampling frequency equal to 8000[samples/sec], the frequency appears to be 1 Khz, while using 16000 sampling frequency the frequency of the signal appear to be 2 Khz.

Which is the correct frequency? And is there other method to calculate the frequency of a signal without a prior knowledge of the sampling frequency?

If you do not have the information about the sampling frequency $F_s$ of your digital data, the best option is to talk about dimensionless relative frequencies $f$, or reduced frequency. The frequencies you observe on periodograms will be $f = F/F_s$, where $F$ would be the true frequency (hoping you have no aliasing). This amounts to saying that your sampling period is $1$ (dimensionless) and that the maximum observable frequency in your signal is $1/2$.

However, it is likely that the actual sampling frequency can be obtained from the data file, the  experiment or sensor, a recorded phenomenon with known frequency (like the 50 or 60 Hz power) or the person who gave you that signal.


What are filter banks What are filter banks and what's their application? What's the meaning of $H(z)$ in descriptions of filter banks?

A filter bank is a collection of bandpass filters designed to split a signal into a number of bands. The centre frequencies of the band pass filters can be spaced linearly or according to a non-linear spacing depending on the intended application.

$H(z)$ denotes the $z$-domain transfer function of a filter in the bank.


DWT versus band-pass filter Which technique gives better result in extracting a band from a signal, let's say the alpha band (8-13 Hz):  using a band-pass-filter or a DWT (Discrete Wavelet Transformation)?

P.S.: I am working in EEG processing.

With high probability, the band-pass filter.

DWT are invertible non-redundant discrete transformations that decompose data onto iterated low- and high-pass filters. Hence, the bands limits are mostly inverses of powers of two. So, you cannot chose arbitrary band values with DWT. More details can be found in the answer to What are the pros and cons of wavelet for filtering compared to conventional filters?

DWT satisfies additional properties (like perfect reconstruction, vanishing moments) that linear filters are not subject to. Hence, the choice in admissible DWT  filters is reduced due to these additional constraints.

Hence, for pure band-pass filtering, a carefully designed pure band-filter can always be better than a DWT (and I usually am pro-wavelet). However, one can design a band-pass filter that is way worse than a DWT processing in certain settings.