Dr Isabelle Niedhammer
email: isabelle.niedhammer@inserm.fr

Simone David

Stéphanie Degioanni

Anne Drummond

Pierre Philip

Occupational Physicians

Drs Acquarone D

F Aicardi

BOUHNIK D, BOULANGER MF, BOULARD J, BOURREAU P, BOURRET D, BOUSTIERE AM, BRETON C, BUGEON G Bordes-Dupuy G, Borrelly N

Canonne Buono-Michel M

CAPELLA D Jf

CERVONI C, CHARRETON D, CHARRIER D, CHAUVIN MA Cavin-Rey M

FABRE D, FAIVRE A, FALLERI R, FERRANDO G Ferrari-Galano J

LOUVAIN S, LOZE M Lotte Mc

Mareel Mc

Workplace bullying and psychotropic drug use: the mediating role of physical and mental health status

Keywords: psychotropic drug use, workplace bullying, mental health, self-reported health OR: odds ratio CI: confidence interval OR adjusted for age, marital status, presence of children, education, occupation, working hours a week, night work, physico-chemical exposures, self-reported health, and depressive symptoms *: p<0.05

Objectives:

The association between workplace bullying and psychotropic drug use is not well established. This study was aimed at exploring the association between workplace bullying, and its characteristics, and psychotropic drug use, and studying the mediating role of physical and mental health. Methods: The study population consisted of a random sample of 3132 men and 4562 women of the working population in the South-East of France.

Workplace bullying, evaluated using the validated instrument elaborated by Leymann, and psychotropic drug use, as well as covariates, were measured using a self-administered questionnaire. Covariates included age, marital status, presence of children, education, occupation, working hours, night work, physico-chemical exposures at work, self-reported health, and depressive symptoms. Statistical analysis was performed using logistic regression analysis, and was carried out separately for men and women. Results: Workplace bullying was strongly associated with psychotropic drug use. Past exposure to bullying increased the risk for this use. The more frequent and the longer the exposure to bullying, the stronger the association with psychotropic drug use. Observing bullying on someone else at the workplace was associated with psychotropic drug use. Adjustment for covariates did not modify the results. Additional adjustment for self-reported health and depressive symptoms reduced the magnitude of the associations, especially for men. Conclusion: The association between bullying and psychotropic drug use was found to be significant and strong, and was partially mediated by physical and mental health.

INTRODUCTION

The use of psychotropic drugs, aimed at bringing a symptomatic response to mental disorders, has increased within the last few decades in industrialised countries, and this category of drugs is now amongst the most commonly prescribed. The ESEMeD project showed that 8% of males and 16% of females took at least one psychotropic drug in the past 12 months in Europe [START_REF] Alonso | Psychotropic drug utilization in Europe: results from the European Study of the Epidemiology of Mental Disorders (ESEMeD) project[END_REF]. In France, 24% of the population (17% of men, 31% of women) were reimbursed for a psychotropic medication during the year 2000 [START_REF] Lecadet | Médicaments psychotropes : consommation et pratiques de prescription en France métropolitaine. I. Données nationales[END_REF] and psychotropic drug use may be higher than in other European countries [START_REF] Alonso | Psychotropic drug utilization in Europe: results from the European Study of the Epidemiology of Mental Disorders (ESEMeD) project[END_REF]. A variety of factors may explain the differences between countries, such as access and utilization of health and mental health services, legal rules for prescription, retail and use of psychotropic drugs, but also actual differences in the prevalence of mental disorders.

Anxiolytics are currently the most commonly prescribed drugs, antidepressants and hypnotics being less frequently used. This pattern is found in Europe in general [START_REF] Alonso | Psychotropic drug utilization in Europe: results from the European Study of the Epidemiology of Mental Disorders (ESEMeD) project[END_REF] and in France in particular [START_REF] Lecadet | Médicaments psychotropes : consommation et pratiques de prescription en France métropolitaine. I. Données nationales[END_REF]. Under-treatment was found to be a common phenomenon: a substantial proportion of subjects with mental disorders did not use psychotropic drugs or were inappropriately treated [START_REF] Alonso | Psychotropic drug utilization in Europe: results from the European Study of the Epidemiology of Mental Disorders (ESEMeD) project[END_REF]. Although benefits are associated with the use of these drugs, concern has also increased regarding their health effects (side effects), especially with long-term use, and the dependence related to use.

Psychotropic drug use may be associated with health, social, and economic consequences, such as cardiac side effects [START_REF] Witchel | Psychotropic drugs, cardiac arrhythmia, and sudden death[END_REF], falls in older age [START_REF] Leipzig | Drugs and falls in older people: a systematic review and meta-analysis: I. Psychotropic drugs[END_REF], and accidents [START_REF] Wadsworth | Psychotropic medication use and accidents, injuries and cognitive failures[END_REF]. Consequently, psychotropic drug use may be a serious public health issue, because of the high prevalence of use, and their consequences, and a better understanding of the factors associated with use may be useful for preventive action.

The causes of psychotropic drug use are complex and certainly multifactorial. Studies reported that psychotropic drug use was associated with female gender, older age, low socioeconomic status, or living alone [START_REF] Alonso | Psychotropic drug utilization in Europe: results from the European Study of the Epidemiology of Mental Disorders (ESEMeD) project[END_REF][START_REF] Lavigne | Psychosocial work environment, interpersonal violence at work and psychotropic drug use among correctional officers[END_REF][START_REF] Lecadet | Médicaments psychotropes : consommation et pratiques de prescription en France métropolitaine. I. Données nationales[END_REF]. To date, occupational factors have been understudied in relation to psychotropic drug use.

Indeed, little research has been undertaken that targets psychotropic drug use in the working population. Although job stress is suspected to be one of the major occupational risk factors for mental health outcomes, the impact of job stress on psychotropic drug use has not been widely studied. Some studies have shown that measures of job stress, such as high psychological demands, low decision latitude, job strain, low social support, iso-strain, low reward, or effort-reward imbalance were associated with psychotropic drug use [START_REF] Lavigne | Psychosocial work environment, interpersonal violence at work and psychotropic drug use among correctional officers[END_REF][START_REF] Moisan | Job strain and psychotropic drug use among white-collar workers[END_REF][START_REF] Pelfrene | Use of benzodiazepine drugs and perceived job stress in a cohort of working men and women in Belgium. Results from the BELSTRESS-study[END_REF][START_REF] Virtanen | Work stress, mental health and antidepressant medication findings from the Health 2000 Study[END_REF]. However, these studies were done on relatively small or selective samples, and/or did not take adequate account of potential confounding factors such as socio-demographic factors, physical and mental health status, or important occupational risk factors such as night work and working hours. Furthermore, a still lower number of studies have explored the association between workplace bullying, considered to be one of the most damaging job stress factors, and psychotropic drug use, and the literature in this area appears still more sparse [START_REF] Appelberg | The use of tranquilizers, hypnotics and analgesics among 18,592 Finnish adults: associations with recent interpersonal conflicts at work or with a spouse[END_REF][START_REF] Lavigne | Psychosocial work environment, interpersonal violence at work and psychotropic drug use among correctional officers[END_REF][START_REF] Richman | Sexual harassment and generalized workplace abuse among university employees: prevalence and mental health correlates[END_REF][START_REF] Vartia | Consequences of workplace bullying with respect to the well-being of its targets and the observers of bullying[END_REF][START_REF] Traweger | Psychosocial factors at the workplace--do they affect substance use? Evidence from the Tyrolean workplace study[END_REF].

Workplace bullying is difficult to evaluate, and no consensus exists regarding its definition.

Here, the definition by Leymann (Leymann, 1996b) was adopted: workplace bullying 'involves hostile and unethical communication, which is directed in a systematic way by one or a few individuals mainly towards one individual who is pushed into a helpless and defenceless position'. Two approaches using self-reported questionnaires have been developed in surveys: inventories of various forms of bullying, and self-report of being exposed to bullying on the basis of a given definition. According to some authors, the combination of both approaches would be adequate to define cases of bullying [START_REF] Einarsen | Bullying at work: epidemiological findings in public and private organizations[END_REF][START_REF] Einarsen | Harassment and bullying at work: a review of the Scandinavian approach[END_REF][START_REF] Vartia | The sources of bullying: psychological work environment and organizational climate[END_REF][START_REF] Vartia | Consequences of workplace bullying with respect to the well-being of its targets and the observers of bullying[END_REF]. Duration and frequency of bullying would also be crucial elements. In the present study, we used the questionnaire elaborated by Leymann, the Leymann Inventory of Psychological Terror (LIPT) (Leymann, 1996a), considered to have the greatest coverage and acceptable reliabilities [START_REF] Cowie | Measuring workplace bullying[END_REF], and also the second approach based on self-report of being exposed to bullying.

The objectives of this study were to examine the associations between workplace bullying and psychotropic drug use. This study attempted to take the previous limitations described earlier into account, as it was based on a large and non-selective sample of the working population, it included a standard measure of exposure to workplace bullying, and detailed information on this exposure, took account of a large number of confounding factors, and studied the specific mediating role of physical and mental health in the association between workplace bullying and psychotropic drug use.

METHODS

Population

This cross-sectional survey was performed by INSERM in 2004 among the general working population in the southeast of France in collaboration with a network of 143 voluntary occupational physicians, who, if working full-time, selected 150 employees each randomly, and invited them to participate in the survey. Occupational medicine is mandatory for all employees in France, consequently, every employee has a medical examination with an occupational physician periodically; at the time of the survey, this was annually. In order to be included in the survey, employees had to have worked for at least 3 months in their company. The survey was based on a self-administered questionnaire which was anonymous, and was returned using a prepaid envelope. Because employees included in the survey were all working at the time of the survey, it could be assumed that those who had a major mental health disorder might be underrepresented in the sample because these people would be more likely to be on sick leave. Several papers have already been published on the topic of workplace bullying using this study sample (Niedhammer et al., 2006b;Niedhammer et al., 2006a;Niedhammer et al., 2007;Niedhammer et al., 2009).

Variables

Our questionnaire included a French version of the LIPT, measuring the experience of 45 forms of bullying within the previous 12 months, as well as frequency and duration of bullying. The reference period for exposure was 12 months making it possible to study exposure to bullying even in a previous job (inclusion criterion in the survey was at least 3 months in the company, consequently exposure to bullying in a previous job was also explored). Afterwards, the employees were given the following definition: 'Bullying may be defined by a situation in which someone is exposed to hostile behaviours on the part of one or more persons in the work environment which aim continually and repeatedly to offend, oppress, maltreat, or to exclude or isolate over a long period of time.' The employees were asked if they perceived themselves as being exposed to bullying within the past 12 months.

Cases of bullying were defined using both Leymann's definition (exposure to at least one form of bullying within the previous 12 months, weekly or more, and for at least 6 months) (Leymann, 1996b) and the self-report of being exposed to bullying, as recommended previously [START_REF] Einarsen | Bullying at work: epidemiological findings in public and private organizations[END_REF][START_REF] Einarsen | Harassment and bullying at work: a review of the Scandinavian approach[END_REF][START_REF] Vartia | The sources of bullying: psychological work environment and organizational climate[END_REF][START_REF] Vartia | Consequences of workplace bullying with respect to the well-being of its targets and the observers of bullying[END_REF]. The psychometric properties of the French version of the LIPT questionnaire were studied in a previous paper (Niedhammer et al., 2006a), and we found that the combined evaluation of bullying increased the convergent and predictive validity compared with Leymann's definition alone. Several other variables were used to characterize the exposure to workplace bullying within the 12 past months, which were the period of exposure (current or past), and the frequency and duration of exposure. In addition, we used a variable describing the fact that the employees may have been observers of bullying experienced by someone else at their workplace within the 12 past months. We also constructed a variable combining the two variables of exposure to bullying and observing bullying by creating four categories: no exposure at all, observer of bullying, exposure to bullying, and both exposure to bullying and observer of bullying.

Psychotropic drug use was measured using one question evaluating whether a medication to sleep, a tranquilizer, or another psychotropic drug for mental disorders had been taken within the previous 12 months. This variable was used as a marker of sleep and mental health disorders. Drug abuse and recreational use of these drugs were not the focus of the study.

Several variables were used as covariates because they are considered as known or suspected risk factors of psychotropic drug use in the literature and may play a confounding role in the association between bullying and psychotropic drug use: age, marital status, presence of children, educational level, occupational groups (derived from the French classification of occupations -INSEE-that is close to the International Standard Classification of Occupations -ISCO-and included: blue collar workers, clerks/service workers, associate professionals, managers/professionals), working hours per week, night work (permanent night work or alternating shifts including night shift), and the number of physico-chemical exposures at work to thermic constraints (outdoor work, cold or hot temperatures), noise, radiation, chemical exposures, or other exposures. Two health-related variables were also studied: poor self-reported health, based on a 4-level scale ranging from 'very good' (coded 1) to 'very poor' (coded 4), and defined by levels 3 and 4, and depressive symptoms measured using the CES-D scale and defined using the available thresholds established for the French population ( 17 for men and  23 for women) to dichotomize the CES-D score [START_REF] Fuhrer | La version française de l'échelle CES-D (Center for Epidemiologic Studies -Depression scale). Description et traduction de l'échelle d'autoévaluation[END_REF]. Age, marital status, and educational level are well-known risk factors of psychotropic drug use [START_REF] Alonso | Psychotropic drug utilization in Europe: results from the European Study of the Epidemiology of Mental Disorders (ESEMeD) project[END_REF][START_REF] Lavigne | Psychosocial work environment, interpersonal violence at work and psychotropic drug use among correctional officers[END_REF][START_REF] Lecadet | Médicaments psychotropes : consommation et pratiques de prescription en France métropolitaine. I. Données nationales[END_REF][START_REF] Moisan | Job strain and psychotropic drug use among white-collar workers[END_REF][START_REF] Pelfrene | Use of benzodiazepine drugs and perceived job stress in a cohort of working men and women in Belgium. Results from the BELSTRESS-study[END_REF]. As our variable of interest was workplace bullying, a careful consideration was also given to potential occupational risk factors of psychotropic drug use, occupational group, working hours, night work, and physico-chemical exposures, that may increase the risk of psychotropic drug use [START_REF] Niedhammer | Psychotropic drug use and shift work among French nurses (1980-1990)[END_REF][START_REF] Pelfrene | Use of benzodiazepine drugs and perceived job stress in a cohort of working men and women in Belgium. Results from the BELSTRESS-study[END_REF].

Finally, self-reported health and depressive symptoms were studied as markers of physical and mental health status, in order to explore their potential mediating role in explaining the association between bullying and psychotropic drug use, something never done before.

Statistical analysis

The crude associations between 6 variables characterizing bullying, i.e. exposure, period, frequency, duration of bullying and the two variables of observing bullying, covariates, and psychotropic drug use were studied using Pearson's Chi-Square test. Logistic regression analysis was then used to adjust for covariates (except self-reported health and depressive symptoms). Six different models were constructed with psychotropic drug use as the dependent variable. In each model, we included as independent variables simultaneously: one of the 6 variables describing bullying, as well as the covariates (model 1). Additional models were also performed with additional adjustment for poor self-reported health and depressive symptoms, that were added in the previous models as covariates (model 2). The mediating role of these two health-related variables was estimated by the change in odds ratio (OR) according to the formula: (OR model 1 -OR model 2 )/(OR model 1 -1) [START_REF] Lynch | Do cardiovascular risk factors explain the relation between socioeconomic status, risk of all-cause mortality, cardiovascular mortality, and acute myocardial infarction?[END_REF]. The OR, or the prevalence-odds ratio as the study design is here cross-sectional, is the ratio of the odds in favour of disease among the exposed to the odds in favour of disease among the unexposed.

As differences may be observed between men and women for the prevalences of occupational exposures (including bullying), and of health outcomes (including psychotropic drug use), and as the associations between exposures and outcomes may also differ between genders, analysis was carried out separately for men and women [START_REF] Niedhammer | How is sex considered in recent epidemiological publications on occupational risks?[END_REF]. Statistical analysis was performed using SAS.

RESULTS

In 2004, 19655 employees were asked to participate in the survey. Among them, 7770 responded to the self-administered questionnaire, leading to a response rate of 40%. Seventy six employees were excluded from the analysis, 57 because they had worked for less than 3 months in their company, and 19 because sex response was missing in the questionnaire.

Thus, the study was based on 7694 employees, 3132 men and 4562 women, with a mean age of 40 (standard deviation: 10.3). The 12-month prevalence of psychotropic drug use was 19% and 33% for men and women respectively.

Leymann's definition alone led to a 12-month prevalence of exposure to bullying of 11% for men and 13% for women. Using the definition of exposure to bullying combining Leymann's definition and the self-reporting of bullying by the employees within the same period, the 12month prevalence of exposure to bullying were respectively 9% and 11% for men and women. These results show that most of those defined as exposed to bullying using Leymann's definition also reported being exposed.

All the associations between the variables of bullying and psychotropic drug use were strongly significant at p<0.001 (Table 1). The prevalence of psychotropic drug use increased among people exposed to workplace bullying, especially among those who were currently exposed. People who were exposed to bullying in the past were also at higher risk of psychotropic drug use than those who had never been exposed. Significant associations were observed between the frequency and duration of exposure to bullying and psychotropic drug use. Observing bullying was also associated with an increase in the prevalence of psychotropic drug use. The highest prevalence of psychotropic drug use was observed for those exposed to bullying (with or without observing it) for both genders.

The prevalence of psychotropic drug use increased with age, among those living alone, among women having no children, among those having a lower educational level, among male clerks or service workers, among night workers, and among women exposed to physicochemical exposures at work. The prevalence of psychotropic drug use also increased strongly with poor self-reported health and depressive symptoms (Table 2).

Table 3 provides the results of logistic regression analysis (model 1). Each model shows the association between each variable of bullying and psychotropic drug use after adjustment for covariates. All these associations were strongly significant at p<0.001, suggesting that covariates did not modify the strong associations observed in Table 1. Exposure to workplace bullying within the last 12 months was found to be a strong risk factor for psychotropic drug use. Past exposure to bullying also increased this risk. The more frequent and the longer the exposure to bullying, the higher the prevalence of use. Tests for trend of frequency and duration of bullying showed dose-response associations significant at p<0.0001. Observing bullying of someone else increased the risk of psychotropic drug use. Exposure to bullying with or without observing bullying at the workplace led to the highest increases in the risk of psychotropic drug use.

These results from logistic regression analysis (Table 3) also showed that older age and living alone for both genders, occupation and night work for men only, and physico-chemical exposures for women only were found to be significant independent risk factors for psychotropic drug use (Table 4).

Additional adjustment for self-reported health and depressive symptoms (model 2) led to a substantial reduction in the magnitude of the ORs, and the associations between bullying and psychotropic drug use remained significant at p<0.001 in women (Table 5). Inclusion of selfreported health and depressive symptoms in the models led to a reduction of the ORs for those exposed to bullying from 3.10 to 1.31 for men and from 3.49 to 1.78 for women, i.e. a reduction of 85% and 69% respectively.

DISCUSSION

Main findings

This study shows that workplace bullying was strongly associated with psychotropic drug use. Dose-response associations were observed between frequency and duration of bullying and psychotropic drug use. Observing bullying of someone else at the workplace was a risk factor of psychotropic drug use. All these associations were independent of potential confounding factors. However, taking physical and mental health into account reduced the magnitude of the associations between bullying and psychotropic drug use, suggesting a mediating role of these health-related variables.

Strengths and limitations of the study

Some limitations of the study deserve to be mentioned. The response rate may be considered as low (40%), but it is similar to previous studies on this sensitive topic [START_REF] Bjorkqvist | Aggression among university employees[END_REF][START_REF] Einarsen | Bullying and harassment at work and their relationships to work environment quality: an exploratory study[END_REF][START_REF] Hoel | The experience of bullying in Great Britain: the impact of organizational status[END_REF][START_REF] Salin | Prevalence and forms of bullying among business professionals: a comparison of two different strategies for measuring bullying[END_REF]. Selection bias may not be ruled out, but the differences between respondents and non-respondents were small for sex, age, economic activities, and occupation. In addition, comparison between the census population and the sample studied suggested that the study population was roughly representative for the same variables (Niedhammer et al., 2006a). A specific type of healthy worker effect (leading to exclude people in poor health from the labour market) may have operated, as this survey did not include employees on sickness absence within the survey period, especially those who were on sick leave because of the health consequences of workplace bullying, leading to a potential underestimation of the association between bullying and psychotropic drug use. The cross-sectional design of our study did not allow us to conclude on the causal nature of the association between workplace bullying and psychotropic drug use, and a reverse causation may not be excluded. A reporting bias may also be suspected as both workplace bullying and psychotropic drug use were measured using self-report. This reporting bias, which is connected to 'common method variance', for example through negative affectivity and social desirability, may lead to inflated associations between the two variables. Another limitation is related to the use of a rather crude measure for psychotropic drug use. Although other authors used a similar measure [START_REF] Lavigne | Psychosocial work environment, interpersonal violence at work and psychotropic drug use among correctional officers[END_REF], this did not allow us to study the specific classes of psychotropics, the reason or duration for medication.

The strengths of this study were the following. Our sample included a very large number of employees of the general working population, allowing us to study a nonselective population as well as men and women separately, which has been shown to be crucial [START_REF] Niedhammer | How is sex considered in recent epidemiological publications on occupational risks?[END_REF]. Another type of healthy worker effect (people in poor health shifted to less exposed jobs) was taken into account as we were able to study exposure to bullying in a previous job.

A validated instrument was used to measure workplace bullying (LIPT), and various variables were constructed to describe the exposure to workplace bullying which has never been done before in the study of psychotropic drug use, and provided detailed information on exposure to bullying and its associations with psychotropic drug use. The prevalence of psychotropic drug use observed for men and women in our study (19-33%) was consistent with previous results in the French population (17-31%) [START_REF] Lecadet | Médicaments psychotropes : consommation et pratiques de prescription en France métropolitaine. I. Données nationales[END_REF], and strong and consistent associations were also found with classical risk factors such as gender, age, and marital status, in accordance with previous studies [START_REF] Alonso | Psychotropic drug utilization in Europe: results from the European Study of the Epidemiology of Mental Disorders (ESEMeD) project[END_REF][START_REF] Lavigne | Psychosocial work environment, interpersonal violence at work and psychotropic drug use among correctional officers[END_REF][START_REF] Lecadet | Médicaments psychotropes : consommation et pratiques de prescription en France métropolitaine. I. Données nationales[END_REF][START_REF] Moisan | Job strain and psychotropic drug use among white-collar workers[END_REF][START_REF] Pelfrene | Use of benzodiazepine drugs and perceived job stress in a cohort of working men and women in Belgium. Results from the BELSTRESS-study[END_REF]. Regarding occupational factors, we found that occupation was significantly associated with psychotropic drug use among men (blue collar workers being at lower risk, and white collar workers at higher risk), in agreement with the results by Pelfrene et al. [START_REF] Pelfrene | Use of benzodiazepine drugs and perceived job stress in a cohort of working men and women in Belgium. Results from the BELSTRESS-study[END_REF]. We observed an association between night work and psychotropic drug use as reported previously [START_REF] Niedhammer | Psychotropic drug use and shift work among French nurses (1980-1990)[END_REF]. An association was also found between physico-chemical exposures at work and psychotropic drug use, that may suggest that other occupational exposures may be risk factors. Finally, as expected, strong and significant associations were found between self-reported health and depressive symptoms, and psychotropic drug use. All these results reinforce the validity of our study. The statistical analysis took into account important covariates in the study of the association between workplace bullying and psychotropic drug use; sociodemographic and occupational factors, as well as health-related factors. Although a substantial number of covariates were taken into account, some variables may have been omitted because they were unavailable in the survey, such as smoking or alcohol consumption. Psychosocial work factors, such as social support, were not taken into account because it is likely that these variables may be risk factors of bullying, and would lead to overadjustment in the analysis. We demonstrated in a previous paper that occupation was associated with bullying, especially that male managers/professionals were less likely to be exposed to bullying (Niedhammer et al., 2007). Our results that were adjusted for occupation (among others) suggested that bullying was associated with psychotropic drug use independently of occupational groups. However, as occupational group may be a risk factor of bullying, adjusting for occupation may also be considered as an overadjustment in the study of the association between bullying and psychotropic drug use. As the association was very strong, adjusting was preferred even if it may be viewed as a conservative strategy.

Whereas sociodemographic and occupational covariates did not modify the strong associations between bullying and psychotropic drug use, health-related covariates reduced these associations, suggesting a mediating role of these covariates.

Comparison with literature

Only very few studies examined the association between workplace bullying and psychotropic drug use. [START_REF] Appelberg | The use of tranquilizers, hypnotics and analgesics among 18,592 Finnish adults: associations with recent interpersonal conflicts at work or with a spouse[END_REF] reported an association between interpersonal conflicts at work and the use of tranquilizers and hypnotics after adjustment for age, social class and neuroticism, and stratified by marital status in Finnish men. [START_REF] Lavigne | Psychosocial work environment, interpersonal violence at work and psychotropic drug use among correctional officers[END_REF] reported an association between psychological harassment and psychotropic drug use after adjustment for gender, age, and several other psychosocial work factors among Canadian correctional officers. The study by [START_REF] Richman | Sexual harassment and generalized workplace abuse among university employees: prevalence and mental health correlates[END_REF] reported an association between generalized workplace abuse and use of prescription drugs (tranquilizers, antidepressants, and sedatives) after adjustment for age, race, and occupation among American male university employees.

The study by [START_REF] Traweger | Psychosocial factors at the workplace--do they affect substance use? Evidence from the Tyrolean workplace study[END_REF] showed that being a victim of bullying was associated with the use of antidepressants, benzodiazepines, analgetics, stimulants, and phytotherapeutics in Austria, but this finding was based on a crude association without adjustment. Vartia [START_REF] Vartia | Consequences of workplace bullying with respect to the well-being of its targets and the observers of bullying[END_REF] reported that bullied employees used both sleep-induced drugs and sedatives more often than the other subjects; observers were also more likely to use them than nonbullied, but the use of these drugs did not correlate with the duration or the frequency of bullying among Finnish municipal employees. However, no information was provided on covariates. In all these studies, the measurement of bullying and/or psychotropic drug use were based on single items. These very few studies also underlined the scarcity of the literature on this topic.

Our results are in agreement with these studies. Furthermore, as never observed before, we found dose-response associations between frequency and duration of bullying and psychotropic drug use. Past exposure to bullying had still an impact on psychotropic drug use, pointing to the long term effects of bullying, even when the exposure had stopped. Being an observer of bullying was a risk factor of psychotropic drug use suggesting that workplaces prone to bullying may have a detrimental effect on employees' health, even if employees are not directly concerned by the phenomenon. Our results are in agreement with the study by Vartia [START_REF] Vartia | Consequences of workplace bullying with respect to the well-being of its targets and the observers of bullying[END_REF] finding an association between observing bullying and the use of sleepinduced drugs and sedatives. Finally, we found that physical and mental health status may mediate the association between workplace bullying and psychotropic drug use, something never reported before. Our study also stated gender differences as physical and mental health status may play a stronger mediating role among men than among women, suggesting that other factors might mediate the association between bullying and psychotropic drug use among women.

Conclusion

Our findings underline the strong association between workplace bullying and psychotropic drug use. Our study also demonstrates that the association between workplace bullying and psychotropic drug use may be mediated by health-related factors, which is a finding that allows a better understanding of the mechanisms leading to psychotropic drug use. More information from prospective etiological studies would be helpful. More research on management styles and organisational policies would also be useful to provide guidelines to policies oriented towards preventing bullying at the workplace. Better understanding of bullying, its determinants and consequences, as well as increasing efforts toward prevention are urgently needed.

Table 2 .

 2 Crude associations between covariates and psychotropic drug use

	Psychotropic drug use
	Men		Women
	Use	Use	Use	Use
	(N)	(%)	(N)	(%)

Table 3 .

 3 Exposure to bullying and psychotropic drug use according to logistic regression

	analysis
			Psychotropic drug use
			Men		Women
		OR	95% CI	OR	95% CI
	Exposure to bullying
	No	1		1
	Yes	3.10	2.36-4.08	3.49	2.86-4.28
	Exposure to bullying
	No exposure	1		1
	Past exposure	2.28	1.10-4.71	3.09	2.14-4.47
	Current exposure	3.22	2.40-4.32	3.56	2.82-4.51
	Frequency of exposure to bullying
	No exposure	1		1
	Weekly	2.30	1.58-3.33	2.96	2.23-3.92
	Daily or almost daily	4.36	2.97-6.40	4.06	3.10-5.33
	Duration of exposure to bullying
	No exposure	1		1
	< 2 y	2.95	1.87-4.64	3.03	2.26-4.06
	 2 y, but < 5 y	2.65	1.75-4.00	3.94	2.85-5.44
	5 years or more	4.25	2.55-7.09	3.84	2.50-5.89
	Observer of bullying
	No	1		1
	Yes	1.92	1.58-2.34	1.92	1.68-2.21
	Combination exposure to bullying/observer
	No exposure	1		1
	Observer	1.71	1.37-2.13	1.55	1.33-1.81
	Exposure to bullying	4.84	2.70-8.66	3.12	2.12-4.61
	Exposure to bullying and observer	3.45	2.52-4.71	4.32	3.41-5.46
	OR: odds ratio CI: confidence interval
	OR adjusted for age, marital status, presence of children, education, occupation, working hours a week, night
	work, and physico-chemical exposures
	Bold OR: significant at 5%
	All bullying variables were significant at p<0.001

Table 4 .

 4 Factors associated with psychotropic drug use according to logistic regression

	analysis
			Psychotropic drug use
			Men		Women
			(N=2911)		(N=4322)
		OR	95% CI	OR	95% CI
	Exposure to bullying	***		***
	No	1		1
	Yes	3.10	2.36-4.08	3.49	2.86-4.28
	Age (years)	***		***
	< 30	1		1
	30-39	1.56	1.12-2.17	1.01	0.81-1.25
	40-49	2.06	1.46-2.90	1.44	1.15-1.78
	50 or more	2.32	1.64-3.29	1.53	1.23-1.90
	Marital status	**		***
	Married, cohabiting	1		1
	Single, separated, divorced, widowed	1.40	1.09-1.78	1.31	1.14-1.51
	Presence of children
	Yes	1		1
	No	1.01	0.80-1.28	1.09	0.93-1.27
	Education
	Primary, lower vocational, lower secondary	1.17	0.89-1.54	1.22	1.01-1.48
	Upper secondary	1.30	0.96-1.77	1.19	0.99-1.43
	University	1		1
	Occupation	*
	Blue collar worker	0.61	0.41-0.92	0.78	0.49-1.24
	Clerks, service workers	1.01	0.70-1.46	1.06	0.78-1.42
	Associate professionals	0.82	0.60-1.11	1.18	0.89-1.56
	Managers, professionals	1		1
	Working hours a week
	<40	1		1
	40	0.94	0.77-1.16	1.08	0.90-1.29
	Night work	*
	No	1		1
	Yes	1.35	1.00-1.83	1.21	0.92-1.60
	Number of physico-chemical exposures			*
	0	1		1
	1	1.20	0.90-1.59	1.06	0.89-1.25
	2	1.49	1.07-2.06	0.98	0.76-1.26
	 3	1.27	0.97-1.67	1.56	1.18-2.08
	OR: odds ratio CI: confidence interval
	*: p<0.05 **: p<0.01 ***: p<0.001
	Bold OR: significant at 5%

Table 5 .

 5 Exposure to bullying and psychotropic drug use according to logistic regression analysis including additional adjustment for self-reported health and depressive symptoms

		Psychotropic drug use
		Men		Women
	OR	95% CI	OR	95% CI

ACKNOWLEDGEMENTS

The authors thank Drs JL Battu, C Beyssier, N Blanchet, C Breton, M Buono-Michel, JF Canonne, G Dalivoust, A Faivre, F Fournier, G Gazazian, G Gibelin-Dol, E Griffaton, T Ibagnes, B Jaubert, J Leroy, M Lozé, G Magallon, P Presseq, G Roux, and C Vitrac, as well as Drs J Chiaroni, C Kaltwasser, and M Signouret, for their help in constructing the network of voluntary occupational physicians, and last but not least S Mocaer and P Sotty for their support and interest in this research program. A pre-test of this survey was performed in 2003 thanks to the contribution of Drs AM Boustière, C Breton, C Cervoni, A Faivre, C Gravier, E Halimi, M Isnard, MC Jacquin, M Lafon-Borelli, J Laudicina, D Londi, M Lozé, G Magallon, V Marcelot, M Méric, C Milliet, F Occhipinti, P Occhipinti, H Raulot-Lapointe, and MO Vincensini. The authors' thanks also go to T Theorell for his help in providing the German version of the LIPT questionnaire, and to I Revue and R Revue for their help in translating and back-translating this questionnaire. Finally, the authors thank all the employees who participated to the survey and made this study possible. This study was supported by the DRTEFP PACA of the French Minister of Labour.


I. INTRODUCTION

E Mpirical Mode Decomposition (EMD) was originally introduced in the late 1990's to study water surface wave evolution [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF]. The EMD can be considered as an emerging technique in signal processing with a very important topic of research and development in various fields such as biomedical signal analysis [START_REF] Hualou | Application of the empirical mode decomposition to the analysis of esophageal manometric data in gastroesophageal reflux disease[END_REF], Hurst exponent estimation [START_REF] Rilling | Empirical mode decomposition, fractional Gaussian noise and Hurst exponent estimation[END_REF], speech processing [START_REF] Chatlani | EMD-based noise estimation and tracking (ENET) with application to speech enhancement[END_REF], texture analysis [START_REF] Nunes | Texture analysis based on local analysis of the bidimensional empirical mode decomposition[END_REF], etc. It decomposes adaptively a given signal, s, into a sum of N AM-FM components, d n (referred to as the Intrinsic Mode Functions, IMFs), plus a residue a N . An IMF is defined [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] as a locally centered function where the number of extrema and the number of zero-crossings must differ at most by one. More precisely, for a given signal s = a 0 , the EMD sequentially computes the N IMFs d n , and N corresponding trends a n , such that a n-1 = a n +d n . The EMD key issue is then the extraction of the N IMFs d n . In practice, such a signal is obtained by stopping a so-called sifting process, using a Cauchy-like criterion [START_REF] Rilling | On empirical mode decomposition and its algorithms[END_REF]. If k denotes the number of iteration in the sifting process, the so-called sifting process can be summarized as follows:

1) Initialization with d n,0 = a n-1 .

. This work was supported by the National Research Agengy (ANR) of France under Grant mv-EMD BLAN07-0314-02.

J. Fleureau, J.-C. Nunes, A. Kachenoura, L. Albera and L. Senhadji are with the INSERM, U642, Rennes F-35000, France, and the Université de Rennes 1, LTSI, Rennes F-35000, France (e-mail: firstname.secondname@univ-rennes1. fr).

2) Computation of the mean envelope M(d n,k ).

3) Extraction of the detail d n,k+1 = d n,k -M(d n,k ). 4) Incrementation of k and return to step 2 if d n,k+1 is not designated as an IMF else stop of the procedure. As proposed by Huang [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF], the mean envelope M(d n,k ) is generally given by the half sum of the upper and the lower envelopes, which are obtained by interpolation between the local maxima points and the local minima points of d n,k , respectively.

Simplicity and efficiency of the original EMD, named Huang hereafter, [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] is seducing but, even if we ignore the theoretical lacks of the classical EMD algorithm, one important limitation lies in its strictly monovariate definition. Indeed, even if the major part of the algorithm seems to be quite general, extension to multivariate signals defined from Ê to Ê D with D > 1 is difficult due to the definition of the mean operator M. Some recent work tried to overcome this limitation. The bivariate approaches developed in [START_REF] Tanaka | Complex empirical mode decomposition[END_REF], [START_REF] Altaf | Rotation invariant complex empirical mode decomposition[END_REF] are extensions of Huang's solution to specifically handle complex-valued signals whereas the algorithm proposed in [START_REF] Rehman | Empirical Mode decomposition for trivariate signals[END_REF] makes use of quaternion to perform trivariate decompositions. Regarding the geometrical algorithms proposed in [START_REF] Rilling | Bivariate empirical mode decomposition[END_REF], they allow us to process bivariate signals. Indeed in [START_REF] Rilling | Bivariate empirical mode decomposition[END_REF], Rilling's algorithm corresponding to the second scheme (called Rilling2 method in the sequel) especially computes a bivariate mean envelope by projecting the signal on different angular planes and by performing a monovariate EMD of each signal projection. Very recently, Rehman et al. [START_REF] Rehman | Multivariate empirical mode decomposition[END_REF] proposed a generalization of Rilling's approach, called Rehman in the sequel, to any multivariate signal: the mean envelope is obtained by averaging multiple signal projections on a regularly sampled hypersphere.

This paper aims at proposing an alternative to Rehman's algorithm [START_REF] Rehman | Multivariate empirical mode decomposition[END_REF] and at unifying mono-and multivariate EMD to process signals with values in Ê D (D ≥ 1) whatever the dimension D is. The proposed method, called 2T-EMD (Turning Tangent EMD), differs by its computational lightness and its algorithmic simplicity. This method is essentially based on a redefinition of the signal mean envelope computed thanks to new characteristic points, which are also defined for D > 1.

The scope of application of the novel algorithm is discussed, and a comparison of the 2T-EMD technique with classical methods is performed on various simulated mono-and multivariate signals. The monovariate behavior of the proposed method on noisy signals is then validated by decomposing fractional Gaussian noises [START_REF] Rilling | Empirical mode decomposition, fractional Gaussian noise and Hurst exponent estimation[END_REF] and an application to surface EEG signal denoising is finally presented.

II. THE 2T-EMD APPROACH FOR BOTH MONO-AND MULTIVARIATE SIGNALS

In order to get a unified framework for mono-and multivariate EMD [START_REF] Fleureau | 3A-EMD: a generalized approach for monovariate and multivariate EMD[END_REF], the signal mean trend and consequently the mean operator M have to be redefined. In some words, the signal mean trend is redefined as the signal which interpolates the barycenters of particular oscillations, called elementary oscillations hereafter.

A. Elementary oscillations and barycenters

An elementary oscillation of a given function s with values in Ê D (D ≥ 1) can be viewed as a piece of s defined between two consecutive local extrema of s. But, for functions with values in a D-dimensional (D > 1) space, the notion of extremum has to be defined.

Let s be a class C 1 function, say differentiable with a continuous first derivative. The function tangent vector to s, denoted by T s , is defined from Ê to Ê D+1 by:

T s : t -→ [1, ds dt (t)] T (1)
Now, let α s be the function given by:

α s : t -→ lim h→0 T s (t -h), T s (t + h) (2)
For every t in Ê, the value α s (t) can be interpreted as the Euclidean inner product of Ê D+1 , denoted by • , • , between the tangents to s just before and after point t. In fact, α s is maximum at point t when both vectors T s (t-h) and T s (t+h) are collinear and it vanishes at point t when both vectors are orthogonal. In other words, function α s is an indicator of the local oscillation of s at every point of Ê. Moreover, due to the continuity of the inner product, we have:

∀t ∈ Ê, α s (t) = lim h→0 T s (t -h), lim h→0 T s (t + h) (3)
Next, since s is in C 1 , we get:

∀t ∈ Ê, α s (t) = T s (t) 2 = 1+ ds dt (t) 2 (4)
where . abusively represents the Euclidean norm of both Ê D+1 and Ê D . Thus, we define an oscillation extremum of function s as a local minimum of function α s . From (4), it also corresponds to a local minimum of the following function:

β s : t -→ β s (t) = ds dt (t) 2 (5)
Clearly, the definition of oscillation extremum does not depend on the output space dimension D, which makes it suitable in multivariate EMD contexts. For D = 1, the reader can check that the oscillation extrema computed from function β s (5) include the classical scalar extrema used in Huang's solution [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF], but also the saddle points (stationary points which are not local extrema) of s and the inflexion points corresponding to positive and negative maxima of the derivative of s. It is noteworthy that oscillation extrema are different of inflexion [START_REF] Delechelle | Empirical mode decomposition: an analytical approach for sifting process[END_REF] and curvature [START_REF] Niang | Décomposition modale empirique: contribution à la modélisation mathématique et application en traitement du signal et de l'image[END_REF] extrema. In addition, curvature extrema may require the computation of the second and third derivatives of s, respectively, which may be more sensitive to sampling

M P1→P2 = t 1 + t 2 2 , 1 t 2 -t 1 t2 t1 s(t) dt T (6)
The mean trend M(s) could then be redefined as the function which interpolates between oscillation barycenters of s. Nevertheless, a straightforward interpolation of all oscillation barycenters of s appears in practice to significantly emphasize the phenomena of over-decomposition of AM-FM signals. A robust computation of the mean trend is preferably obtained for 2T-EMD by averaging two envelopes: a first envelope interpolates the even indexed barycenters which include signal borders, and a second envelope interpolates the odd indexed barycenters which also include signal borders. This operation lightly increases the algorithm computational complexity (see subsection II-D for more details) but makes its behavior considerably more robust.

Indeed, this robustness may be roughly justified in a monodimensional context. The spline interpolation typically behaves as a linear filter (up to some border effects) whose frequency response I(f ), for unit spaced knots, slowly decreases in the interval f ∈ [0, 1] and approximately vanishes when f > 1 [START_REF] Unser | Splines: A perfect fit for signal processing[END_REF]. Let's consider now a narrowband AM-FM signal whose carrier frequency is f c . If the modulation is not too strong, the typical spacing between oscillation barycenters for such a signal is 1/(2f c ). In the case of a direct interpolation between oscillation barycenters, the spline interpolation could then have a non-negligible content for f ∈ [0, 2f c ], which includes f c and therefore allows the interpolation to oscillate at the same frequency as the AM-FM signal. On the other hand, the use of two interleaved spline interpolations leads to a frequency content typically limited to f < f c , which approximately guarantees that the resulting mean trend oscillates more slowly than the initial AM-FM signal. Therefore, the use of two envelopes prevents from over-decomposition phenomena and also has the advantage to bring the frequency resolution of the proposed method on a par with the traditional EMD approaches (for signals with simple waveforms). It finally approximately guarantees that the local frequency of successive IMFs decreases. Those remarks seem to be also empirically verified (see section 2) in a multivariate context.

Consequently, the resulting mean operator M enjoys a frequency property similar to that of the classical one but can be applied to both mono-and multivariate signals without any dimension restriction.

B. Implementation of the 2T-EMD technique

Several important points have to be considered to achieve a robust and efficient implementation of the 2T-EMD algorithm. First, the derivative necessary to the computation of M(s) is calculated by means of a centered finite difference scheme. However, in practice the resulting numerical derivative may have more local minima than the theoretical one. Minima are thus detected using an ε precision equal to 10 -15 in our code. More precisely, the i-th time index t i will be a minimum of

β s if β s (t i ) + ε ≤ β s (t i-1 ) and β s (t i ) + ε ≤ β s (t i+1
). Therefore, this ε precision limits the influence of numerical artifacts for the extrema detection but one should be careful and adjust the ε value in the case of very high sampling frequency to avoid the miss of legitimate extrema. The integral involved in the computation of barycenters is then calculated using a simple rectangle method. Second, the interpolation is performed using cubic splines with classical boundary conditions as for the classical EMD, where the signal borders are directly added to the list of estimated oscillation barycenters. In other words, there is no complex management of boundary conditions. Third, the sifting process is stopped using a modified Cauchylike criterion. More especially, the following normalized Eu-

clidean norm ||d n,k+1 (t) -d n,k (t)||/||d n,k (t)|| is computed at each point t.
If a given percentage of the latter norms, for instance 90%, is below a certain threshold fixed to 10 -2 in our code, then the sifting process is stopped. The normalized Euclidean norm could be obviously not defined for d n,k (t) = 0 but the set of points where this zero value is reached is likely to have a zero Lebesgue measure which justifies in practice the use of such a criterion.

C. Scope of application

Let's now specify in more details the set of signals that 2T-EMD can decompose successfully. As previously mentioned, the considered signals have to be in class C 1 or at least, in the case of irregular signals (especially real and/or noisy signals), an appropriate numerical estimation of the derivative has then to be proposed (see section II-B). Note that some existing methods based on inflexion [START_REF] Delechelle | Empirical mode decomposition: an analytical approach for sifting process[END_REF] or curvature [START_REF] Niang | Décomposition modale empirique: contribution à la modélisation mathématique et application en traitement du signal et de l'image[END_REF] extrema require more regularity. In addition, if {d n } 1≤n≤N represents the theoretical set of IMFs composing the signal to analyze, it is obvious that any d n with a piecewise constant function β dn [START_REF] Nunes | Texture analysis based on local analysis of the bidimensional empirical mode decomposition[END_REF] is not visible by 2T-EMD (assuming that the definition of derivative is extended to irregular signal as mentioned before). Indeed a piecewise constant function has no ε-local minimum. Thus, any signal in class C 1 having one IMF with a piecewise constant derivative norm cannot be properly decomposed by 2T-EMD. It mainly concerns monovariate signals with piecewise linear IMFs, and bivariate signals with purely circular rotating IMFs. Nevertheless, any bivariate rotating signal with a sufficient eccentricity can be decomposed by 2T-EMD.

D. Note on computational complexity

The computational complexity of the 2T-EMD algorithm can be precisely evaluated and compared with the one of classical methods. For a given EMD algorithm, let N , K n , d n,k be the number of extracted IMFs, the number of sifting iterations performed to extract the n-th IMF and the n-th IMF computed at the k-th iteration of the sifting process, respectively. In addition, M H (n, k), M R (n, k, p), M RM (n, k, p), and M 2T (d n,k ) will denote the number of extrema detected in d n,k by Huang [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF], the number of extrema detected in the p-th projection of d n,k by Rilling2 [START_REF] Rilling | Bivariate empirical mode decomposition[END_REF] when P projection planes are used, the number of extrema detected in the pth projection of d n,k by Rehman [START_REF] Rehman | Multivariate empirical mode decomposition[END_REF] when P projection directions are used, and the number of barycenters detected in d n,k by 2T-EMD. Then the number F (d n,k+1 ) of multiplications and divisions (usually called number of flops) of one sifting iteration necessary to obtain d n,k+1 from d n,k is given in table I for the four methods. These results were obtained by considering a standard tridiagonal implementation of the spline interpolation and a signal from Ê to Ê D of L samples. The total computational cost, C(s), of the four methods can be obtained straightforwardly by summing the elementary complexities, F (d n,k ), given in table I over both the number of iterations and the number of IMFs.

Method D Numerical complexity F (d n,k+1 ) Huang [1] 1 18L + 15M H (n, k) Rilling2 [10] 2 L(11P + 2) + 15 P /2 p=1 M R (n, k, p) Rehman [11] 2 LP (2D + 18) + 15 P p=1 M RM (n, k, p) 2T-EMD AE * D(19L + 16M 2T (n, k)) + M 2T (n, k)

III. SIMULATION RESULTS

The aim of this section is to analyze the performance of 2T-EMD on several classes of simulated mono-and multivariate signals satisfying the conditions given in subsection II-C. First the stability and the convergence speed of 2T-EMD sifting process are studied and compared with three classical methods: i) Huang [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] for monovariate signals (D = 1), ii) Rilling2 [START_REF] Rilling | Bivariate empirical mode decomposition[END_REF] using P = 8 projection planes for bivariate signals (D = 2) and iii) Rehman [START_REF] Rehman | Multivariate empirical mode decomposition[END_REF] for trivariate and quadrivariate signals using P = 2D projection directions on the associated hypersphere. From an implementation point of view, the sifting process termination criterion and the border management of those classical techniques are identical to the ones used in 2T-EMD (see subsection II-B for more details). Secondly, a benchmark involving those classical EMD algorithms is presented to evaluate and compare the 2T-EMD performance in the context of a full decomposition. Finally, an example of quadrivariate decomposition illustrates the ability of 2T-EMD to work for a signal dimension greater than 3. First, let's introduce the signal selection and the performance criteria used in this section.

A. Signals selection

The is given by:

d (12) 1 (t) = 0.3 exp(0.23(1 + t)) sin(100πt) d (12) 2
(t) = exp(0.23(1 + t)) sin(50πt)

The third signal s 13 = d

(13) 1 + d (13) 2
+ a (13) (a (13) represents a residue) is defined by:

d (13) 1 (t) = sin(200π(1+t) 2 )/2 d (13) 2 (t) = 2(1 + t) sin(100πt + 1.5) a (13) (t) = 5t 2
Eventually, the fourth signal s 14 = d

(14) 1 + d (14) 2 + d (14) 3
with a residue a (14) is given by:

d (14) 1 (t) = 2 sin(250πt) d (14) 2 (t) = 3 cos(π(1.7t + 7.3) 2 ) d (14) 3 (t) = exp(0.23(1 + t)) cos(π(2.58t + 21.95) 2 ) a (14) (t) = 3t
More particularly, signal s 11 is the sum of two sinusoidal components. Next s 12 is the sum of two sinusoidal components modulated in amplitude (from 1 to 2 on interval T ) with extrema of oscillation that do not match with classical extrema. As far as s 13 is concerned, it is the sum of one is given by: d is defined by: used to show the efficiency of the approach in higher dimensions is given by: d

d (22) 1 (t) = [exp(t) cos(π/4) sin(80πt + 1.5), exp (

B. Performance criteria

The different EMD algorithms analyzed in this section have been evaluated and compared in terms of performance and numerical complexity. Given I a subinterval of T = [-1, 2], let's define the following quadratic errors:

e I (d n ) = I ||d n (t) -dn (t)|| 2 dt I ||d n (t)|| 2 dt , e I (s) = N n=1 e I (d n ) (7)
where dn denotes the estimate of the n-th IMF d n of signal s.

The latter errors allow us to evaluate the ability of the EMD algorithms to accurately extract one or all expected IMFs. By considering I = [0, 1] where border effects should be low, both errors allow us to compute a performance independent of any border effects. On the contrary, with I = [-1, 0] ∪ [1, 2], both errors permit to evaluate the ability of a given algorithm to minimize border effects (all algorithms have the same border management in our study). Indeed, such effects are often critical especially in real life data and their minimization should facilitate the practical exploitation of the IMFs. Eventually, the numerical complexity of an EMD method is evaluated using the criterion C(s) presented in section II-D.

C. Sifting process analysis

In this first experiment, 2T-EMD's sifting process is studied and compared to that of Huang, Rilling2 and Rehman in terms of convergence. To do so, the sifting process is launched for 100 iterations to extract the first IMF of all monovariate and multivariate signals presented in section III-A. Figure 1 displays criterion e [0,1] (d 1 ) at the output of 2T-EMD, Huang, Rilling2 and Rehman as a function of the number of iterations of the sifting process. On each curve, a circular-shaped marker indicates the iteration for which the Cauchy-like sifting stop criterion would have in practice interrupted the sifting process. First, one can observe that the sifting process of 2T-EMD converges for almost all monovariate and bivariate test signals with precision and speed very similar to those obtained by the considered classical algorithms. The proposed algorithm even succeeds in extracting the first component of signal s 12 whereas the standard Huang algorithm failed in such a decomposition. Except for the multivariate signals s 24 and s 41 , this simulation tends to show the stability of the proposed algorithm during a sifting process. The signal s 24 with four IMFs including especially two FM components seems to suffer from over-decomposition. It is noteworthy that this sensitivity is also observed in Rilling2 and may especially suggest that the first component d is not considered as an IMF by both algorithms. Regarding the signal s 41 , Rehman seems to be less sensitive to over-decomposition phenomenon than 2T-EMD. However, for all proposed signals, including signals s 24 and s 41 , the sifting stop criterion prevents from any overdecomposition phenomenon by stopping the sifting process at an appropriate optimal iteration. In addition, note that the quadratic error e [0,1] (d 1 ) at marker points is satisfactory for all studied signals (the maximum error encountered in this simulation is equal to 0.02) and is, globally, comparable to the one obtained by the classical algorithms.

D. Performance study of the full process

This subsection aims at comparing the full process of 2T-EMD for D ∈ {1, 2, 3, 4} with Huang, Rilling2 and Rehman using criteria e [0,1] (s) and e [-1,0]∪ [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF][START_REF] Hualou | Application of the empirical mode decomposition to the analysis of esophageal manometric data in gastroesophageal reflux disease[END_REF] (s) and C(s). Results are depicted in figure 2: the left, middle and right columns deal with the monovariate, bivariate and multivariate (D ∈ {3, 4}) cases, respectively. It appears that for no border effects both performance and computational complexity of 2T-EMD are globally equivalent to those of the reference methods. More finely and as specified previously, in the monovariate case, signal s 12 with extrema of oscillation which differ from the standard extrema, is better decomposed by our approach. On the contrary, the signal s 13 with a first fast FM IMF seems to be better processed by Huang. In the bivariate case, Rilling2 seems to offer a better performance on the nearly pure rotating signal s 21 but 2T-EMD is more efficient on signal s 22 made of rotated planar components. For signals s 31 and s 41 , the behaviors of 2T-EMD and Rehman are quasi-similar. Regarding border effects on performance, 2T-EMD generally offers slightly more efficient results and seems to provide a more accurate management of border effects than reference approaches. When focusing on the computational complexity, figure 2 clearly shows that for D = 1 the computational complexities of 2T-EMD and Huang are quasi-equivalent. On the contrary, for D = 2, 2T-EMD generally requires less sifting iterations and less computational operations than Rilling2. For D > 2, the results clearly show the computational efficiency of 2T-EMD, which is about three times cheaper than Rehman. From a more illustrative point of view, figure 3 (a) represents the quadrivariate signal s 41 (dark line) and the associated local mean (gray line) on the restricted [0.50; 0.52] time interval projected on the three frames (W, X, Y ), (X, Y, Z) and (Y, Z, W ). The local mean seems to nicely go through the original signal on the three frames. Figure 3 (b) presents the expected and computed decompositions of signal s 41 projected on the four main axis (from left to right), namely W, X, Y Note that only a central zoom of the IMFs on the temporal axis has been represented for the sake of clarity. A good behavior of the algorithm can be observed and the low residue highlights the efficiency of the proposed method. This result shows clearly that, for signal s 41 , the sifting stop criterion prevents from any over-decomposition phenomenon.

In conclusion, all these results exhibit the ability of 2T-EMD to process signals with values in Ê D for any dimension D ≥ 1. They also show the advantage of the use of oscillation barycenters in a general multivariate context.

IV. DECOMPOSITION OF FRACTIONAL GAUSSIAN NOISE

One important property, enjoyed by the classical Huang algorithm, is the filter bank property observed when decomposing a fractional Gaussian noise (fGn). This section aims at demonstrating how the 2T-EMD algorithm also enjoys this filter bank property by reproducing simulations and results already obtained in previous works with the standard Huang's algorithm [START_REF] Wu | A study of the characteristics of white noise using the empirical mode decomposition method[END_REF], [START_REF] Rilling | Empirical mode decomposition, fractional Gaussian noise and Hurst exponent estimation[END_REF]. Note that computing 2T-EMD on such a noisy signal may be feasible using an adapted definition of the derivative (see II-B). As in [START_REF] Rilling | Empirical mode decomposition, fractional Gaussian noise and Hurst exponent estimation[END_REF], extensive simulations are carried out on fGn's, with Hurst coefficient H varying from 0.1 to 0.9. The data length is typically set to 512 and for each value of H, 1000 independent sample paths of fGn are generated via the algorithm described in [START_REF] Perrin | Synthèse exacte et efficace du mouvement brownien fractionnaire 1D[END_REF] 2] and suggest the hierarchical structure of an equivalent filter bank as shown in [START_REF] Rilling | Empirical mode decomposition, fractional Gaussian noise and Hurst exponent estimation[END_REF] and [START_REF] Wu | A study of the characteristics of white noise using the empirical mode decomposition method[END_REF] for the classical EMD. For all IMFs d n (with n > 1), a self-similarity in this filter-bank could also be further checked showing that 2T-EMD approximately acts on fGn as a dyadic filter bank of constant-Q bandpass filters for high values of H (H ≥ 0.5). Assuming this filter bank structure, and as shown in [START_REF] Rilling | Empirical mode decomposition, fractional Gaussian noise and Hurst exponent estimation[END_REF] for the classical EMD, it becomes possible to get access to the Hurst exponent via the variance progression across IMFs by var(d n ) ∝ ρ V. APPLICATION TO SURFACE EEG DENOISING Surface ElectroEncephaloGraphy (EEG) is a popular neuroimaging technique used for exploring human brain activity. While this technique is simple and low cost, the obtained signals suffer from noise and artifacts, such as broken wire contacts, ocular movements (ElectroOculoGram, EOG), muscular activity, etc. Thus, one of the challenging tasks in signal processing is to detect and extract very weak non-stationary brain source signals corrupted by noise and artifacts from EEG data. Such issues are of great interest when EEG features are used for diagnosis and assistance [START_REF] Kachenoura | ICA : a potential tool for BCI systems[END_REF]. Some recent studies show that the use of classical EMD [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] in order to denoise EEG data [START_REF] Rutkowski | Clustering of spectral patterns based on EMD components of EEG channels with applications to neurophysiological signals separation[END_REF], [START_REF] Molla | Separation of EOG artifacts from EEG signals using bivariate EMD[END_REF] and/or to detect some EEG patterns [START_REF] Rutkowski | EMD approach to multichannel EEG data -the amplitude and phase components clustering analysis[END_REF], give interesting results. In this section, we propose to remove EOG artifacts from contaminated EEG data by using the proposed 2T-EMD algorithm. Note that, the goal here is not to propose a new method in order to denoise the surface EEG data, but to show the behavior of 2T-EMD, in comparison to that of the existing EMD method, in the case of real world signals. The used EEG signals are issued from our polysomnographic database [START_REF] Poree | Blind source separation for ambulatory sleep recording[END_REF]. More precisely, EEG signals (figure 5, line 1) are acquired from two temporal electrodes, in front of the higher part of the ears, denoted by F7m and F8m (where "m" stands for modified, see [START_REF] Poree | Blind source separation for ambulatory sleep recording[END_REF] for more details). Additionally, EOG reference recordings (figure 5, line 6) are taken from two temporal sensors located near each eye (EOGL and EOGR), slightly moved toward the median plan in order to simultaneously observe horizontal and vertical eye movements.

Rilling2 and the bivariate 2T-EMD method are then applied with a limited number of 10 IMFs, for each bivariate observation [F7m,F8m] and [EOGL,EOGR], respectively. The obtained bivariate results are then projected on two main axis (left head side and right head side), namely F7p and F8p for EEG data and EOGLp and EOGRp in the case of EOG data ("p" stands for projected). In order to identify the projected IMFs related to the EOG artifacts in the bivariate EEG observation [F7m,F8m], we firstly compute the Fourier transform of all projected IMFs extracted from the two bivariate observations ([F7m,F8m] and [EOGR,EOGL]) [START_REF] Rutkowski | Clustering of spectral patterns based on EMD components of EEG channels with applications to neurophysiological signals separation[END_REF]. Then, for each bivariate method and each brain hemisphere (F7p/EOGLp on the left side and F8p/EOGRp on the right side), i) oneminus-correlation distance is computed for between projected IMFs in order to evaluate any spectral similarity between the IMFs issued from EEG signals [F7m,F8m] and EOG reference signals [EOGR,EOGL], and ii) the distances are hierarchically clustered using the single linkage algorithm [START_REF] Rutkowski | EMD approach to multichannel EEG data -the amplitude and phase components clustering analysis[END_REF]. It is noteworthy that, the clustering procedure first treats the power spectrum vector of each projected IMF as a singleton cluster and then successively aggregates both most similar clusters, until all clusters merge into a single cluster that contains the power spectrum vectors of all projected IMFs. For each method and for each projected plan, the set of IMFs for which the distances (one-minus-correlation) below 0.2 are considered as similar and originate from EOG artifacts. The sum of those IMFs [F7p EOG 2T-EMD , F8p EOG 2T-EMD , F7p EOG Rilling2 , F8p EOG Rilling2 ] are then compared in figure 5 lation coefficients greater than 0.92). This high correlation is also observed between the residual high-frequency channels (associated to EEG activity in our case) obtained by 2T-EMD method, [F7p EEG 2T-EMD , F8p EEG 2T-EMD ], and those recovered by Rilling2 namely, [F7p EEG Rilling2 , F8p EEG Rilling2 ]. All these results tend to show that 2T-EMD has similar behavior in comparison to classical EMD algorithms [START_REF] Rilling | Empirical mode decomposition, fractional Gaussian noise and Hurst exponent estimation[END_REF] in the case of real world data. In addition, they demonstrate that the 2T-EMD technique which takes advantage of the mutual-information contained in the two signals, preserves the frequency information of each channel.

VI. CONCLUSION AND PERSPECTIVES

A novel approach for EMD computation is proposed in this paper. The algorithm called 2T-EMD is mainly based on a new geometric definition of the mean envelope operator. Under certain assumptions on signals, this novel definition enables both mono-and multivariate decompositions without any modification in the 2T-EMD algorithm, without any signal projection and with light computational cost. These two last points are the key benefits regarding the existing approaches of the literature. Details to obtain a robust implementation of 2T-EMD have been listed and justified. In addition, simulations and comparisons performed in this work suggest that 2T-EMD seems to offer a satisfactory convergence during the sifting step and a good robustness to over-decomposition. Competitive performance and computational complexities plus a good border management also characterize 2T-EMD. Simulations on noisy signals and an application to real data show that the proposed algorithm has quasi-identical behavior in comparison with classical EMD algorithms and demonstrate the interest of 2T-EMD in practical situations. Conditions that signals should verify to be successfully decomposed by 2T-EMD have been precisely enumerated in order to help the user. It would be now more particularly interesting to evaluate the performance of 2T-EMD in other real life contexts for which classical EMD approaches have already demonstrated interesting properties and to compare it with that of other techniques [START_REF] Kachenoura | ICA : a potential tool for BCI systems[END_REF].

Turning

  Tangent Empirical Mode Decomposition: a framework for mono-and multivariate signals Julien Fleureau, Jean-Claude Nunes, Amar Kachenoura, Laurent Albera, Member, IEEE, and Lotfi Senhadji, Senior Member, IEEE Abstract-A novel Empirical Mode Decomposition (EMD) algorithm, called 2T-EMD, for both mono-and multivariate signals is proposed in this paper. It differs from the other approaches by its computational lightness and its algorithmic simplicity. The method is essentially based on a redefinition of the signal mean envelope, computed thanks to new characteristic points, which offers the possibility to decompose multivariate signals without any projection. The scope of application of the novel algorithm is specified, and a comparison of the 2T-EMD technique with classical methods is performed on various simulated mono-and multivariate signals. The monovariate behaviour of the proposed method on noisy signals is then validated by decomposing a fractional Gaussian noise and an application to real life EEG data is finally presented. Index Terms-Mono-and Multivariate Empirical Mode Decomposition, Intrinsic Mode functions, Analysis of non-linear and non-stationary signals, Hurst exponent estimation, Extrema and barycenters of oscillation, Filter bank structure, EEG denoising, Time varying representation.

  test signals are mono-, bi-, tri-and quadrivariate signals with various AM-FM behaviors, defined on the time interval T = [-1, 2] and sampled at f s = 10 kHz. Monovariate signals. Four monovariate signals of the form s 1i = n d (1i) n are used in this study where d (1i) n is the n-th IMF of s 1i . The first signal s 11 = d

  FM component (from 200 to 400 Hz on interval [0, 1]), one linear AM component and a quadratic residue. Eventually s 14 is the sum of one sinusoidal component, one FM component (from 10 to 20 Hz on interval T ), one AM (from 1 to 2 on interval T ) -FM (from 50 to 70 Hz on interval T ) component and a linear residue. Bivariate signals. Four bivariate signals of the form s 2i = n d (2i) n are used in this analysis where d (2i) n is the n-th IMF of s 2i . The first signal s 21 = d

  = [1.4 cos(40πt), 1.7 sin(40πt)] T The second signal s 22 = d

  = [3 sin(500πt + 2), 3.5 sin(500πt), 2 cos(500πt), exp(0.23t -0.46) sin(500πt)] T d (41) 2 (t) = [0.5 cos(120πt + 1.2), 5 exp (0.23(1 + t)) cos(π(3.16t + 25.3) 2 ), 2 exp(0.23(1 + t)) sin(π(3.16t + 25.3) 2 ), 0.5 sin(130πt + 1.2)] T d AM-FM IMFs of signal s 41 .

Fig. 1 .

 1 Fig. 1. 2T-EMD quadratic error on first IMF as a function of the number of iterations during the sifting process.

Fig. 2 .

 2 Fig. 2. Comparative study of 2T-EMD versus Huang (1D), Rilling2 (2D) and Rehman (3D, 4D) reference methods.

Fig. 4 .

 4 Fig.4. 2T-EMD applied to fGn decompositions. At the left, IMF zerocrossings: the (base 2) logarithm of the average number of zero-crossings is plotted as a function of the IMF number, for 3 different values of the Hurst exponent (H: 0.2 (circles), 0.5 (crosses), and 0.8 (stars)). At the right, IMF variance and estimation of the Hurst exponent H: the IMF log-variance slope plotted as a function of the Hurst exponent H and the associated linear regression over 0.5 ≤ H ≤ 0.9.

  (α-1)n H with the specific choice α = 2H -1. When plotted as a function of the Hurst exponent H, the IMF log-variance slope is almost linear when H ≥ 0.5, in accordance with the simplified model p(H) = 2log(ρ H )(H -1) as depicted in the right column of figure 4. Those results consequently highlight the very similar properties shared by 2T-EMD and Huang ([3, figure 4]), and may validate the use of the proposed algorithm on irregular signals.

Fig. 3 .

 3 Fig. 3. 2T-EMD quadrivariate decomposition of the signal s 41 : (a) original signal (dark) and its 2T-EMD local mean (gray), and (b) exact (dashed line) and estimated (plain line) IMFs projected on the four main axis W, X, Y and Z.

Fig. 5 .

 5 Fig. 5. Comparison between the bivariate 2T-EMD approach and Rilling2 in the context of EOG artifacts removing from EEG.

  In fact, signals s 21 , s 23 and s 24 are globally rotating signals with various AM and FM modulations. Regarding the conditions addressed in subsection II-C the latter signals have a large enough eccentricity to be processed by 2T-EMD. As far as signal s 22 is concerned, it is the result of planar components after a rotation around the temporal axis.

					(24) 1	+ d (24) 2	+ d (24) 3	is given by:
		d (24) 1	(t) = [0.3 cos(400πt), 0.8 sin(500πt)] T
		d (24) 2	(t) = [4 cos(100πt), 7 sin(100πt)] T
		d (24) 3	(t) = [5 exp(0.23t -0.46) cos(π(2.24t + 4.47) 2 ),
					2 exp(0.23t -0.46) sin(π(2.24t + 4.47) 2 )] T
		Trivariate signal. One trivariate signal s 31 = d (31) 1 d (31) 3 is used for comparison and is defined by:	+ d (31) 2	+
		d (31) 1	(t) = [sin(540πt), 2 sin(560πt), 1.5 sin(540πt)] T
		d (31) 2	(t) = [exp (0.14(1 + t)) cos(200πt), 2 cos(200πt),
					2 sin(200πt + 1.2)] T
		d (31) 3	(t) = 3 exp (0.16t + 1.07) cos(π(1.83t + 7.30) 2 ),
					2 cos(π(1.83t + 7.30) 2 ),
					4 exp (0.16t + 1.07) sin(π(1.83t + 7.30) 2 ) T
		where d (31) 1	, d (31) 2	and d (31) 3	denote the three AM-FM IMFs
		of signal s 31 .
		Quadrivariate signal. The quadrivariate signal s 41 = d (41) 1 d (41) 2 (41) + d 3	+
		t) sin(π/4) sin(80πt + 1.5)] T
	d (22) 2	(t) = [cos(π/4) sin(2π(2.5(1 + t)) (23) 1 + d (23) 2 is given by:
	(23) d 1

2

), sin(π/4) sin(2π(2.5(1 + t)) 2

)] T The third signal s 23 = d (t) = [exp(0.23t -0.46)(0.3 cos(π(2.58t + 14.2) 2 ), 0.9 sin(π(2.58t + 14.2) 2 ))] T d (23) 2 (t) = [0.4 sin(20πt), 0.7 cos(24πt)] T The fourth signal s 24 = d

  and decomposed using 2T-EMD on a limited number of 7 IMFs (it should be mentioned that the sifting process is stopped when 90% of values ||d n,k+1 (t) -d n,k (t)||/||d n,k (t)|| are lower than 10 -2 , see section II-B). For each value of H, the number of zero-crossings z H [n] in the n-th mode (1 ≤ n ≤ 7) is firstly evaluated. A linear regression of the mean log number of zero-crossings log(z H [n]) on the mode number n is then computed. The good fitness of such a regression, represented in the left column of figure 4 for each value of H, strongly suggests that the number of zero-crossings z H [n], which is a rough indication of the mean frequency of each mode n, is a decreasing exponential function of the mode number, i.e., z H [n] ∝ ρ -n H with ρ H close to 2 . These results are very similar to those obtained in [3, figure

VII. ACKNOWLEDGMENT

We would like to thank all anonymous reviewers for their valuable and constructive comments for improving the quality of the paper. This work was supported by the French Government under an ANR contract, namely mv-EMD (BLAN07-0314-02).